!> \file gfdl_cloud_microphys.F90
!! This file contains the full GFDL cloud microphysics (Chen and Lin (2013)
!! \cite chen_and_lin_2013 and Zhou et al. 2019 \cite zhou_etal_2019).
!! The module is paired with 'gfdl_fv_sat_adj', which performs the "fast"
!! processes
!>author Shian-Jiann Lin, Linjiong Zhou
!***********************************************************************
!*                   GNU Lesser General Public License
!*
!* This file is part of the GFDL Cloud Microphysics.
!*
!* The GFDL Cloud Microphysics is free software: you can
!* redistribute it and/or modify it under the terms of the
!* GNU Lesser General Public License as published by the
!* Free Software Foundation, either version 3 of the License, or
!* (at your option) any later version.
!*
!* The GFDL Cloud Microphysics is distributed in the hope it will be
!* useful, but WITHOUT ANYWARRANTY; without even the implied warranty
!* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
!* See the GNU General Public License for more details.
!*
!* You should have received a copy of the GNU Lesser General Public
!* License along with the GFDL Cloud Microphysics.
!* If not, see <http://www.gnu.org/licenses/>.
!***********************************************************************
! =======================================================================
!>\defgroup mod_gfdl_cloud_mp GFDL Cloud MP modules
!!\ingroup gfdlmp
!! This module contains the column GFDL Cloud microphysics scheme.
module gfdl_cloud_microphys_mod

    ! use mpp_mod, only: stdlog, mpp_pe, mpp_root_pe, mpp_clock_id, &
    ! mpp_clock_begin, mpp_clock_end, clock_routine, &
    ! input_nml_file
    ! use diag_manager_mod, only: register_diag_field, send_data
    ! use time_manager_mod, only: time_type, get_time
    ! use constants_mod, only: grav, rdgas, rvgas, cp_air, hlv, hlf, pi => pi_8
    ! use fms_mod, only: write_version_number, open_namelist_file, &
    ! check_nml_error, file_exist, close_file

   use module_mp_radar

   implicit none

   private

   public gfdl_cloud_microphys_mod_driver, gfdl_cloud_microphys_mod_init, &
          gfdl_cloud_microphys_mod_end, cloud_diagnosis
!   public wqs1, wqs2, qs_blend, wqsat_moist, wqsat2_moist
!   public qsmith_init, qsmith, es2_table1d, es3_table1d, esw_table1d
!   public setup_con, wet_bulb

   real :: missing_value = - 1.e10

   logical :: module_is_initialized = .false.
   logical :: qsmith_tables_initialized = .false.

   character (len = 17) :: mod_name = 'gfdl_cloud_microphys'

   real, parameter :: n0r = 8.0e6, n0s = 3.0e6, n0g = 4.0e6
   real, parameter :: rhos = 0.1e3, rhog = 0.4e3
   real,                 parameter :: grav = 9.80665       !< gfs: acceleration due to gravity
   real,                 parameter :: rdgas = 287.05       !< gfs: gas constant for dry air
   real,                 parameter :: rvgas = 461.50       !< gfs: gas constant for water vapor
   real,                 parameter :: cp_air = 1004.6      !< gfs: heat capacity of dry air at constant pressure
   real,                 parameter :: hlv = 2.5e6          !< gfs: latent heat of evaporation
   real,                 parameter :: hlf = 3.3358e5       !< gfs: latent heat of fusion
   real,                 parameter :: pi = 3.1415926535897931 !< gfs: ratio of circle circumference to diameter

   ! real, parameter :: rdgas = 287.04                     !< gfdl: gas constant for dry air

   ! real, parameter :: cp_air = rdgas * 7. / 2.           ! 1004.675, heat capacity of dry air at constant pressure
   real, parameter :: cp_vap = 4.0 * rvgas                 !< 1846.0, heat capacity of water vapore at constnat pressure
   ! real, parameter :: cv_air = 717.56                    ! satoh value
   real, parameter :: cv_air = cp_air - rdgas              !< 717.55, heat capacity of dry air at constant volume
   ! real, parameter :: cv_vap = 1410.0                    ! emanuel value
   real, parameter :: cv_vap = 3.0 * rvgas                 !< 1384.5, heat capacity of water vapor at constant volume

   ! the following two are from emanuel's book "atmospheric convection"
   ! real, parameter :: c_ice = 2106.0                     ! heat capacity of ice at 0 deg c: c = c_ice + 7.3 * (t - tice)
   ! real, parameter :: c_liq = 4190.0                     ! heat capacity of water at 0 deg c

   real, parameter :: c_ice = 1972.0                       !< gfdl: heat capacity of ice at - 15 deg c
   real, parameter :: c_liq = 4185.5                       !< gfdl: heat capacity of water at 15 deg c
   ! real, parameter :: c_liq = 4218.0                     ! ifs: heat capacity of liquid at 0 deg c

   real, parameter :: eps = rdgas / rvgas                  !< 0.6219934995
   real, parameter :: zvir = rvgas / rdgas - 1.            !< 0.6077338443

   real, parameter :: t_ice = 273.16                       !< freezing temperature
   real, parameter :: table_ice = 273.16                   !< freezing point for qs table

   ! real, parameter :: e00 = 610.71                       ! gfdl: saturation vapor pressure at 0 deg c
   real, parameter :: e00 = 611.21                         !< ifs: saturation vapor pressure at 0 deg c

   real, parameter :: dc_vap = cp_vap - c_liq              !< - 2339.5, isobaric heating / cooling
   real, parameter :: dc_ice = c_liq - c_ice               !< 2213.5, isobaric heating / colling

   real, parameter :: hlv0 = hlv                           !< gfs: evaporation latent heat coefficient at 0 deg c
   ! real, parameter :: hlv0 = 2.501e6                     ! emanuel appendix - 2
   real, parameter :: hlf0 = hlf                           !< gfs: fussion latent heat coefficient at 0 deg c
   ! real, parameter :: hlf0 = 3.337e5                     ! emanuel

   real, parameter :: lv0 = hlv0 - dc_vap * t_ice          !< 3.13905782e6, evaporation latent heat coefficient at 0 deg k
   real, parameter :: li00 = hlf0 - dc_ice * t_ice         !< - 2.7105966e5, fusion latent heat coefficient at 0 deg k

   real, parameter :: d2ice = dc_vap + dc_ice              !< - 126, isobaric heating/cooling
   real, parameter :: li2 = lv0 + li00                     !< 2.86799816e6, sublimation latent heat coefficient at 0 deg k

   real, parameter :: qrmin = 1.e-8                        !< min value for rain
   real, parameter :: qvmin = 1.e-20                       !< min value for water vapor (treated as zero)
   real, parameter :: qcmin = 1.e-12                       !< min value for cloud condensates

   real, parameter :: vr_min = 1.e-3                       !< min fall speed for rain
   real, parameter :: vf_min = 1.e-5                       !< min fall speed for cloud ice, snow, graupel

   real, parameter :: dz_min = 1.e-2                       !< use for correcting flipped height

   real, parameter :: sfcrho = 1.2                         !< surface air density
   real, parameter :: rhor = 1.e3                          !< density of rain water, lin83
    ! intercept parameters

   real, parameter :: rnzr = 8.0e6 ! lin83
   real, parameter :: rnzs = 3.0e6 ! lin83
   real, parameter :: rnzg = 4.0e6 ! rh84
   real, parameter :: rnzh = 4.0e4 ! lin83 --- lmh 29 sep 17

    ! density parameters

   real, parameter :: rhoh = 0.917e3 ! lin83 --- lmh 29 sep 17

   public rhor, rhos, rhog, rhoh, rnzr, rnzs, rnzg, rnzh
   real :: cracs, csacr, cgacr, cgacs, csacw, craci, csaci, cgacw, cgaci, cracw !< constants for accretions
   real :: acco (3, 4)                                     !< constants for accretions
   real :: cssub (5), cgsub (5), crevp (5), cgfr (2), csmlt (5), cgmlt (5)

   real :: es0, ces0
   real :: pie, rgrav, fac_rc
   real :: c_air, c_vap

   real :: lati, latv, lats, lat2, lcp, icp, tcp           !< used in Bigg mechanism and wet bulk

   real :: d0_vap                                          !< the same as dc_vap, except that cp_vap can be cp_vap or cv_vap
   real :: lv00                                            !< the same as lv0, except that cp_vap can be cp_vap or cv_vap

   ! cloud microphysics switchers

   integer :: icloud_f = 0                                 !< cloud scheme
   integer :: irain_f = 0                                  !< cloud water to rain auto conversion scheme

   logical :: de_ice = .false.                             !< to prevent excessive build - up of cloud ice from external sources
   logical :: sedi_transport = .true.                      !< transport of momentum in sedimentation
   logical :: do_sedi_w = .false.                          !< transport of vertical motion in sedimentation
   logical :: do_sedi_heat = .true.                        !< transport of heat in sedimentation
   logical :: prog_ccn = .false.                           !< do prognostic ccn (yi ming's method)
   logical :: do_qa = .true.                               !< do inline cloud fraction
   logical :: rad_snow = .true.                            !< consider snow in cloud fraciton calculation
   logical :: rad_graupel = .true.                         !< consider graupel in cloud fraction calculation
   logical :: rad_rain = .true.                            !< consider rain in cloud fraction calculation
   logical :: fix_negative = .false.                       !< fix negative water species
    logical :: do_setup = .true. !< setup constants and parameters
   logical :: p_nonhydro = .false.                         !< perform hydrosatic adjustment on air density

   real, allocatable :: table (:), table2 (:), table3 (:), tablew (:)
   real, allocatable :: des (:), des2 (:), des3 (:), desw (:)

    logical :: tables_are_initialized = .false.

   ! logical :: master
   ! integer :: id_rh, id_vtr, id_vts, id_vtg, id_vti, id_rain, id_snow, id_graupel, &
   ! id_ice, id_prec, id_cond, id_var, id_droplets
   real, parameter :: dt_fr = 8.                           !< homogeneous freezing of all cloud water at t_wfr - dt_fr
   ! minimum temperature water can exist (moore & molinero nov. 2011, nature)
   ! dt_fr can be considered as the error bar

   real :: p_min = 100.                                    !< minimum pressure (pascal) for mp to operate

   ! slj, the following parameters are for cloud - resolving resolution: 1 - 5 km

   ! qi0_crt = 0.8e-4
   ! qs0_crt = 0.6e-3
   ! c_psaci = 0.1
   ! c_pgacs = 0.1

   ! -----------------------------------------------------------------------
   ! namelist parameters
   ! -----------------------------------------------------------------------

   real :: cld_min = 0.05                                  !< minimum cloud fraction
   real :: tice = 273.16                                   !< set tice = 165. to trun off ice - phase phys (kessler emulator)

   real :: t_min = 178.                                    !< min temp to freeze - dry all water vapor
   real :: t_sub = 184.                                    !< min temp for sublimation of cloud ice
   real :: mp_time = 150.                                  !< maximum micro - physics time step (sec)

   ! relative humidity increment

   real :: rh_inc = 0.25                                   !< rh increment for complete evaporation of cloud water and cloud ice
   real :: rh_inr = 0.25                                   !< rh increment for minimum evaporation of rain
   real :: rh_ins = 0.25                                   !< rh increment for sublimation of snow

   ! conversion time scale

   real :: tau_r2g = 900.                                  !< rain freezing during fast_sat
   real :: tau_smlt = 900.                                 !< snow melting
   real :: tau_g2r = 600.                                  !< graupel melting to rain
   real :: tau_imlt = 600.                                 !< cloud ice melting
   real :: tau_i2s = 1000.                                 !< cloud ice to snow auto-conversion
   real :: tau_l2r = 900.                                  !< cloud water to rain auto-conversion
   real :: tau_v2l = 150.                                  !< water vapor to cloud water (condensation)
   real :: tau_l2v = 300.                                  !< cloud water to water vapor (evaporation)
   real :: tau_g2v = 900.                                  !< graupel sublimation
   real :: tau_v2g = 21600.                                !< graupel deposition -- make it a slow process

   ! horizontal subgrid variability

   real :: dw_land = 0.20                                  !< base value for subgrid deviation / variability over land
   real :: dw_ocean = 0.10                                 !< base value for ocean

   ! prescribed ccn

   real :: ccn_o = 90.                                     !< ccn over ocean (cm^ - 3)
   real :: ccn_l = 270.                                    !< ccn over land (cm^ - 3)

   real :: rthresh = 10.0e-6                               !< critical cloud drop radius (micro m)

   ! -----------------------------------------------------------------------
   ! wrf / wsm6 scheme: qi_gen = 4.92e-11 * (1.e3 * exp (0.1 * tmp)) ** 1.33
   ! optimized: qi_gen = 4.92e-11 * exp (1.33 * log (1.e3 * exp (0.1 * tmp)))
   ! qi_gen ~ 4.808e-7 at 0 c; 1.818e-6 at - 10 c, 9.82679e-5 at - 40c
   ! the following value is constructed such that qc_crt = 0 at zero c and @ - 10c matches
   ! wrf / wsm6 ice initiation scheme; qi_crt = qi_gen * min (qi_lim, 0.1 * tmp) / den
   ! -----------------------------------------------------------------------

   real :: sat_adj0 = 0.90                                 !< adjustment factor (0: no, 1: full) during fast_sat_adj

   real :: qc_crt = 5.0e-8                                 !< mini condensate mixing ratio to allow partial cloudiness

   real :: qi_lim = 1.                                     !< cloud ice limiter to prevent large ice build up

   real :: ql_mlt = 2.0e-3                                 !< max value of cloud water allowed from melted cloud ice
   real :: qs_mlt = 1.0e-6                                 !< max cloud water due to snow melt

   real :: ql_gen = 1.0e-3                                 !< max cloud water generation during remapping step if fast_sat_adj = .t.
   real :: qi_gen = 1.82e-6                                !< max cloud ice generation during remapping step

   ! cloud condensate upper bounds: "safety valves" for ql & qi

   real :: ql0_max = 2.0e-3                                !< max cloud water value (auto converted to rain)
   real :: qi0_max = 1.0e-4                                !< max cloud ice value (by other sources)

   real :: qi0_crt = 1.0e-4                                !< cloud ice to snow autoconversion threshold (was 1.e-4);
                                                           !! qi0_crt is highly dependent on horizontal resolution
   real :: qr0_crt = 1.0e-4                                !< rain to snow or graupel/hail threshold
                                                           ! lfo used * mixing ratio * = 1.e-4 (hail in lfo)
   real :: qs0_crt = 1.0e-3                                !< snow to graupel density threshold (0.6e-3 in purdue lin scheme)

   real :: c_paut = 0.55                                   !< autoconversion cloud water to rain (use 0.5 to reduce autoconversion)
   real :: c_psaci = 0.02                                  !< accretion: cloud ice to snow (was 0.1 in zetac)
   real :: c_piacr = 5.0                                   !< accretion: rain to ice:
   real :: c_cracw = 0.9                                   !< rain accretion efficiency
   real :: c_pgacs = 2.0e-3                                !< snow to graupel "accretion" eff. (was 0.1 in zetac)

   ! decreasing clin to reduce csacw (so as to reduce cloud water --- > snow)

   real :: alin = 842.0                                    !< "a" in lin1983
   real :: clin = 4.8                                      !< "c" in lin 1983, 4.8 -- > 6. (to ehance ql -- > qs)

   ! fall velocity tuning constants:

   logical :: const_vi = .false.                           !< if .t. the constants are specified by v * _fac
   logical :: const_vs = .false.                           !< if .t. the constants are specified by v * _fac
   logical :: const_vg = .false.                           !< if .t. the constants are specified by v * _fac
   logical :: const_vr = .false.                           !< if .t. the constants are specified by v * _fac

   ! good values:

   real :: vi_fac = 1.                                     !< if const_vi: 1 / 3
   real :: vs_fac = 1.                                     !< if const_vs: 1.
   real :: vg_fac = 1.                                     !< if const_vg: 2.
   real :: vr_fac = 1.                                     !< if const_vr: 4.

   ! upper bounds of fall speed (with variable speed option)

   real :: vi_max = 0.5                                    !< max fall speed for ice
   real :: vs_max = 5.0                                    !< max fall speed for snow
   real :: vg_max = 8.0                                    !< max fall speed for graupel
   real :: vr_max = 12.                                    !< max fall speed for rain

   ! cloud microphysics switchers

   logical :: fast_sat_adj = .false.                       !< has fast saturation adjustments
   logical :: z_slope_liq = .true.                         !< use linear mono slope for autocconversions
   logical :: z_slope_ice = .false.                        !< use linear mono slope for autocconversions
   logical :: use_ccn = .false.                            !< must be true when prog_ccn is false
   logical :: use_ppm = .false.                            !< use ppm fall scheme
   logical :: mono_prof = .true.                           !< perform terminal fall with mono ppm scheme
   logical :: mp_print = .false.                           !< cloud microphysics debugging printout
   logical :: do_hail = .false.                            !< use hail parameters instead of graupel

   ! real :: global_area = - 1.

   real :: log_10, tice0, t_wfr

    integer :: reiflag = 1
    ! 1: Heymsfield and Mcfarquhar, 1996
    ! 2: Wyser, 1998
     
    logical :: tintqs = .false. !< use temperature in the saturation mixing in PDF 

    real :: rewmin = 5.0, rewmax = 10.0
    real :: reimin = 10.0, reimax = 150.0
    real :: rermin = 10.0, rermax = 10000.0
    real :: resmin = 150.0, resmax = 10000.0
    real :: regmin = 300.0, regmax = 10000.0
    
   ! -----------------------------------------------------------------------
   ! namelist
   ! -----------------------------------------------------------------------

   namelist / gfdl_cloud_microphysics_nml /                                  &
       mp_time, t_min, t_sub, tau_r2g, tau_smlt, tau_g2r, dw_land, dw_ocean, &
       vi_fac, vr_fac, vs_fac, vg_fac, ql_mlt, do_qa, fix_negative, vi_max,  &
       vs_max, vg_max, vr_max, qs_mlt, qs0_crt, qi_gen, ql0_max, qi0_max,    &
       qi0_crt, qr0_crt, fast_sat_adj, rh_inc, rh_ins, rh_inr, const_vi,     &
       const_vs, const_vg, const_vr, use_ccn, rthresh, ccn_l, ccn_o, qc_crt, &
       tau_g2v, tau_v2g, sat_adj0, c_piacr, tau_imlt, tau_v2l, tau_l2v,      &
       tau_i2s, tau_l2r, qi_lim, ql_gen, c_paut, c_psaci, c_pgacs,           &
       z_slope_liq, z_slope_ice, prog_ccn, c_cracw, alin, clin, tice,        &
       rad_snow, rad_graupel, rad_rain, cld_min, use_ppm, mono_prof,         &
       do_sedi_heat, sedi_transport, do_sedi_w, de_ice, icloud_f, irain_f,   &
       mp_print, reiflag, rewmin, rewmax, reimin, reimax, rermin, rermax,    &
       resmin, resmax, regmin, regmax, tintqs, do_hail

   public                                                                    &
       mp_time, t_min, t_sub, tau_r2g, tau_smlt, tau_g2r, dw_land, dw_ocean, &
       vi_fac, vr_fac, vs_fac, vg_fac, ql_mlt, do_qa, fix_negative, vi_max,  &
       vs_max, vg_max, vr_max, qs_mlt, qs0_crt, qi_gen, ql0_max, qi0_max,    &
       qi0_crt, qr0_crt, fast_sat_adj, rh_inc, rh_ins, rh_inr, const_vi,     &
       const_vs, const_vg, const_vr, use_ccn, rthresh, ccn_l, ccn_o, qc_crt, &
       tau_g2v, tau_v2g, sat_adj0, c_piacr, tau_imlt, tau_v2l, tau_l2v,      &
       tau_i2s, tau_l2r, qi_lim, ql_gen, c_paut, c_psaci, c_pgacs,           &
       z_slope_liq, z_slope_ice, prog_ccn, c_cracw, alin, clin, tice,        &
       rad_snow, rad_graupel, rad_rain, cld_min, use_ppm, mono_prof,         &
       do_sedi_heat, sedi_transport, do_sedi_w, de_ice, icloud_f, irain_f,   &
       mp_print, reiflag, rewmin, rewmax, reimin, reimax, rermin, rermax,    &
       resmin, resmax, regmin, regmax, tintqs, do_hail 

contains

! -----------------------------------------------------------------------
! the driver of the gfdl cloud microphysics
! -----------------------------------------------------------------------

!>\ingroup mod_gfdl_cloud_mp
!! This subroutine is the driver of the GFDL cloud microphysics
subroutine gfdl_cloud_microphys_mod_driver (                                    &
            iis, iie, jjs, jje, kks, kke, ktop, kbot,                           &
            qv, ql, qr, qi, qs, qg, qa, qn,                                     &
            qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt, qa_dt, pt_dt, pt, w,      &
            uin, vin, udt, vdt, dz, delp, area, dt_in, land,                    &
            rain, snow, ice, graupel, hydrostatic, phys_hydrostatic,            &
            p, lradar, refl_10cm, reset, pfils, pflls)

   implicit none

   ! Interface variables
   integer, intent (in) :: iis, iie, jjs, jje ! physics window
   integer, intent (in) :: kks, kke ! vertical dimension
   integer, intent (in) :: ktop, kbot ! vertical compute domain

   real, intent (in) :: dt_in ! physics time step

   real, intent (in), dimension (iis:iie, jjs:jje) :: area ! cell area
   real, intent (in), dimension (iis:iie, jjs:jje) :: land ! land fraction

   real, intent (in), dimension (iis:iie, jjs:jje, kks:kke) :: delp, dz, uin, vin
   real, intent (in), dimension (iis:iie, jjs:jje, kks:kke) :: pt, qv, ql, qr, qg, qa, qn

   real, intent (inout), dimension (iis:iie, jjs:jje, kks:kke) :: qi, qs
   real, intent (inout), dimension (iis:iie, jjs:jje, kks:kke) :: pt_dt, qa_dt, udt, vdt, w
   real, intent (inout), dimension (iis:iie, jjs:jje, kks:kke) :: qv_dt, ql_dt, qr_dt
   real, intent (inout), dimension (iis:iie, jjs:jje, kks:kke) :: qi_dt, qs_dt, qg_dt

   real, intent (out), dimension (iis:iie, jjs:jje) :: rain, snow, ice, graupel

   logical, intent (in) :: hydrostatic, phys_hydrostatic

   !integer, intent (in) :: seconds
   real, intent (in), dimension (iis:iie, jjs:jje, kks:kke) :: p
   logical, intent (in) :: lradar
   real, intent (out), dimension (iis:iie, jjs:jje, kks:kke) :: refl_10cm
   logical, intent (in) :: reset
   real, intent (out), dimension (iis:iie, jjs:jje, kks:kke) :: pfils, pflls

   ! Local variables
   logical :: melti = .false.

   real :: mpdt, rdt, dts, convt, tot_prec

   integer :: i, j, k
   integer :: is, ie, js, je ! physics window
   integer :: ks, ke ! vertical dimension
   integer :: days, ntimes, kflip

   real, dimension (iie-iis+1, jje-jjs+1) :: prec_mp, prec1, cond, w_var, rh0

   real, dimension (iie-iis+1, jje-jjs+1,kke-kks+1) :: vt_r, vt_s, vt_g, vt_i, qn2

   real, dimension (size(pt,1), size(pt,3)) :: m2_rain, m2_sol

   real :: allmax
!+---+-----------------------------------------------------------------+
!For 3D reflectivity calculations
   real, dimension(ktop:kbot):: qv1d, t1d, p1d, qr1d, qs1d, qg1d, dBZ
!+---+-----------------------------------------------------------------+

   is = 1
   js = 1
   ks = 1
   ie = iie - iis + 1
   je = jje - jjs + 1
   ke = kke - kks + 1
   ! call mpp_clock_begin (gfdl_mp_clock)

   ! -----------------------------------------------------------------------
   ! define heat capacity of dry air and water vapor based on hydrostatical property
   ! -----------------------------------------------------------------------

   if (phys_hydrostatic .or. hydrostatic) then
       c_air = cp_air
       c_vap = cp_vap
       p_nonhydro = .false.
   else
       c_air = cv_air
       c_vap = cv_vap
       p_nonhydro = .true.
   endif
   d0_vap = c_vap - c_liq
   lv00 = hlv0 - d0_vap * t_ice

   if (hydrostatic) do_sedi_w = .false.

   ! -----------------------------------------------------------------------
   ! define latent heat coefficient used in wet bulb and Bigg mechanism
   ! -----------------------------------------------------------------------

   latv = hlv
   lati = hlf
   lats = latv + lati
   lat2 = lats * lats

   lcp = latv / cp_air
   icp = lati / cp_air
   tcp = (latv + lati) / cp_air

   ! tendency zero out for am moist processes should be done outside the driver

   ! -----------------------------------------------------------------------
   ! define cloud microphysics sub time step
   ! -----------------------------------------------------------------------

   mpdt   = min (dt_in, mp_time)
   rdt    = 1. / dt_in
   ntimes = nint (dt_in / mpdt)

   ! small time step:
   dts = dt_in / real (ntimes)

   ! call get_time (time, seconds, days)

   ! -----------------------------------------------------------------------
   ! initialize precipitation
   ! -----------------------------------------------------------------------

   do j = js, je
       do i = is, ie
           graupel (i, j) = 0.
           rain (i, j) = 0.
           snow (i, j) = 0.
           ice (i, j) = 0.
           cond (i, j) = 0.
       enddo
   enddo

   pfils = 0.
   pflls = 0.

   ! -----------------------------------------------------------------------
   ! major cloud microphysics
   ! -----------------------------------------------------------------------

   do j = js, je
       call mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, qg,&
           qa, qn, dz, is, ie, js, je, ks, ke, ktop, kbot, j, dt_in, ntimes,  &
           rain (:, j), snow (:, j), graupel (:, j), ice (:, j), m2_rain,     &
           m2_sol, cond (:, j), area (:, j), land (:, j), udt, vdt, pt_dt,    &
           qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt, qa_dt, w_var, vt_r,      &
           vt_s, vt_g, vt_i, qn2)
       do k = ktop, kbot
           do i = is, ie
               pfils(i, j, k) = m2_sol (i, k)
               pflls(i, j, k) = m2_rain(i, k)
           enddo
       enddo
   enddo

   ! -----------------------------------------------------------------------
   ! no clouds allowed above ktop
   ! -----------------------------------------------------------------------

   if (ks < ktop) then
       do k = ks, ktop
           if (do_qa) then
               do j = js, je
                   do i = is, ie
                       qa_dt (i, j, k) = 0.
                   enddo
               enddo
           else
               do j = js, je
                   do i = is, ie
                       ! qa_dt (i, j, k) = - qa (i, j, k) * rdt
                       qa_dt (i, j, k) = 0. ! gfs
                   enddo
               enddo
           endif
       enddo
   endif

   ! -----------------------------------------------------------------------
   ! diagnostic output
   ! -----------------------------------------------------------------------

   ! if (id_vtr > 0) then
   ! used = send_data (id_vtr, vt_r, time, is_in = iis, js_in = jjs)
   ! endif

   ! if (id_vts > 0) then
   ! used = send_data (id_vts, vt_s, time, is_in = iis, js_in = jjs)
   ! endif

   ! if (id_vtg > 0) then
   ! used = send_data (id_vtg, vt_g, time, is_in = iis, js_in = jjs)
   ! endif

   ! if (id_vti > 0) then
   ! used = send_data (id_vti, vt_i, time, is_in = iis, js_in = jjs)
   ! endif

   ! if (id_droplets > 0) then
   ! used = send_data (id_droplets, qn2, time, is_in = iis, js_in = jjs)
   ! endif

   ! if (id_var > 0) then
   ! used = send_data (id_var, w_var, time, is_in = iis, js_in = jjs)
   ! endif

   ! convert to mm / day

   convt = 86400. * rdt * rgrav
   do j = js, je
       do i = is, ie
           rain (i, j) = rain (i, j) * convt
           snow (i, j) = snow (i, j) * convt
           ice (i, j) = ice (i, j) * convt
           graupel (i, j) = graupel (i, j) * convt
           prec_mp (i, j) = rain (i, j) + snow (i, j) + ice (i, j) + graupel (i, j)
       enddo
   enddo

   ! if (id_cond > 0) then
   ! do j = js, je
   ! do i = is, ie
   ! cond (i, j) = cond (i, j) * rgrav
   ! enddo
   ! enddo
   ! used = send_data (id_cond, cond, time, is_in = iis, js_in = jjs)
   ! endif

   ! if (id_snow > 0) then
   ! used = send_data (id_snow, snow, time, iis, jjs)
   ! used = send_data (id_snow, snow, time, is_in = iis, js_in = jjs)
   ! if (mp_print .and. seconds == 0) then
   ! tot_prec = g_sum (snow, is, ie, js, je, area, 1)
   ! if (master) write (*, *) 'mean snow = ', tot_prec
   ! endif
   ! endif
   !
   ! if (id_graupel > 0) then
   ! used = send_data (id_graupel, graupel, time, iis, jjs)
   ! used = send_data (id_graupel, graupel, time, is_in = iis, js_in = jjs)
   ! if (mp_print .and. seconds == 0) then
   ! tot_prec = g_sum (graupel, is, ie, js, je, area, 1)
   ! if (master) write (*, *) 'mean graupel = ', tot_prec
   ! endif
   ! endif
   !
   ! if (id_ice > 0) then
   ! used = send_data (id_ice, ice, time, iis, jjs)
   ! used = send_data (id_ice, ice, time, is_in = iis, js_in = jjs)
   ! if (mp_print .and. seconds == 0) then
   ! tot_prec = g_sum (ice, is, ie, js, je, area, 1)
   ! if (master) write (*, *) 'mean ice_mp = ', tot_prec
   ! endif
   ! endif
   !
   ! if (id_rain > 0) then
   ! used = send_data (id_rain, rain, time, iis, jjs)
   ! used = send_data (id_rain, rain, time, is_in = iis, js_in = jjs)
   ! if (mp_print .and. seconds == 0) then
   ! tot_prec = g_sum (rain, is, ie, js, je, area, 1)
   ! if (master) write (*, *) 'mean rain = ', tot_prec
   ! endif
   ! endif
   !
   ! if (id_rh > 0) then !not used?
   ! used = send_data (id_rh, rh0, time, iis, jjs)
   ! used = send_data (id_rh, rh0, time, is_in = iis, js_in = jjs)
   ! endif
   !
   !
   ! if (id_prec > 0) then
   ! used = send_data (id_prec, prec_mp, time, iis, jjs)
   ! used = send_data (id_prec, prec_mp, time, is_in = iis, js_in = jjs)
   ! endif

   ! if (mp_print) then
   ! prec1 (:, :) = prec1 (:, :) + prec_mp (:, :)
   ! if (seconds == 0) then
   ! prec1 (:, :) = prec1 (:, :) * dt_in / 86400.
   ! tot_prec = g_sum (prec1, is, ie, js, je, area, 1)
   ! if (master) write (*, *) 'daily prec_mp = ', tot_prec
   ! prec1 (:, :) = 0.
   ! endif
   ! endif

   ! call mpp_clock_end (gfdl_mp_clock)
    if(lradar) then
       ! Only set melti to true at the output times
       if (reset) then
         melti=.true.
       else
         melti=.false.
       endif
       do j = js, je
          do i = is, ie
             do k = ktop,kbot
               kflip = kbot-ktop+1-k+1
               t1d(k)  = pt(i,j,kflip)
               p1d(k)  = p(i,j,kflip)
               qv1d(k) = qv(i,j,kflip)/(1-qv(i,j,kflip))
               qr1d(k) = qr(i,j,kflip)
               qs1d(k) = qs(i,j,kflip)
               qg1d(k) = qg(i,j,kflip)
             enddo
             call refl10cm_gfdl (qv1d, qr1d, qs1d, qg1d,                 &
                       t1d, p1d, dBZ, ktop, kbot, i, j, melti)
             do k = ktop,kbot
                kflip = kbot-ktop+1-k+1
                refl_10cm(i,j,kflip) = MAX(-35., dBZ(k))
             enddo
          enddo
       enddo
    endif

end subroutine gfdl_cloud_microphys_mod_driver

! -----------------------------------------------------------------------
!>\ingroup mod_gfdl_cloud_mp
!>\brief GFDL cloud microphysics, major program, and is based on
!!  Lin et al.(1983) \cite lin_et_al_1983 and
!! Rutledge and Hobbs (1984) \cite rutledge_and_hobbs_1984.
!!
!>\section detmpdrv GFDL Cloud mpdrv General Algorithm
subroutine mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs,     &
        qg, qa, qn, dz, is, ie, js, je, ks, ke, ktop, kbot, j, dt_in, ntimes, &
        rain, snow, graupel, ice, m2_rain, m2_sol, cond, area1, land,         &
        u_dt, v_dt, pt_dt, qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt, qa_dt,   &
        w_var, vt_r, vt_s, vt_g, vt_i, qn2)

    implicit none

    logical, intent (in) :: hydrostatic

    integer, intent (in) :: j, is, ie, js, je, ks, ke
    integer, intent (in) :: ntimes, ktop, kbot

    real, intent (in) :: dt_in

    real, intent (in), dimension (is:) :: area1, land

    real, intent (in), dimension (is:, js:, ks:) :: uin, vin, delp, pt, dz
    real, intent (in), dimension (is:, js:, ks:) :: qv, ql, qr, qg, qa, qn

    real, intent (inout), dimension (is:, js:, ks:) :: qi, qs
    real, intent (inout), dimension (is:, js:, ks:) :: u_dt, v_dt, w, pt_dt, qa_dt
    real, intent (inout), dimension (is:, js:, ks:) :: qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt

    real, intent (inout), dimension (is:) :: rain, snow, ice, graupel, cond

    real, intent (out), dimension (is:, js:) :: w_var

    real, intent (out), dimension (is:, js:, ks:) :: vt_r, vt_s, vt_g, vt_i, qn2

    real, intent (out), dimension (is:, ks:) :: m2_rain, m2_sol

    real, dimension (ktop:kbot) :: qvz, qlz, qrz, qiz, qsz, qgz, qaz
    real, dimension (ktop:kbot) :: vtiz, vtsz, vtgz, vtrz
    real, dimension (ktop:kbot) :: dp0, dp1, dz0, dz1
    real, dimension (ktop:kbot) :: qv0, ql0, qr0, qi0, qs0, qg0, qa0
    real, dimension (ktop:kbot) :: t0, den, den0, tz, p1, denfac
    real, dimension (ktop:kbot) :: ccn, c_praut, m1_rain, m1_sol, m1
    real, dimension (ktop:kbot) :: u0, v0, u1, v1, w1

    real :: cpaut, rh_adj, rh_rain
    real :: r1, s1, i1, g1, rdt, ccn0
    real :: dt_rain, dts
    real :: s_leng, t_land, t_ocean, h_var
    real :: cvm, tmp, omq
    real :: dqi, qio, qin

    integer :: i, k, n

    dts = dt_in / real (ntimes)
    dt_rain = dts * 0.5
    rdt = 1. / dt_in

    ! -----------------------------------------------------------------------
    ! use local variables
    ! -----------------------------------------------------------------------
    
    !GJF: assign values to intent(out) variables that are commented out
    w_var = 0.0
    vt_r = 0.0
    vt_s = 0.0
    vt_g = 0.0
    vt_i = 0.0
    qn2 = 0.0
    !GJF
    
    do i = is, ie

        do k = ktop, kbot
            qiz (k) = qi (i, j, k)
            qsz (k) = qs (i, j, k)
        enddo

        ! -----------------------------------------------------------------------
        !> - Prevent excessive build-up of cloud ice from external sources.
        ! -----------------------------------------------------------------------

        if (de_ice) then
            do k = ktop, kbot
                qio = qiz (k) - dt_in * qi_dt (i, j, k) ! original qi before phys
                qin = max (qio, qi0_max) ! adjusted value
                if (qiz (k) > qin) then
                    qsz (k) = qsz (k) + qiz (k) - qin
                    qiz (k) = qin
                    dqi = (qin - qio) * rdt ! modified qi tendency
                    qs_dt (i, j, k) = qs_dt (i, j, k) + qi_dt (i, j, k) - dqi
                    qi_dt (i, j, k) = dqi
                    qi (i, j, k) = qiz (k)
                    qs (i, j, k) = qsz (k)
                endif
            enddo
        endif

        do k = ktop, kbot

            t0 (k) = pt (i, j, k)
            tz (k) = t0 (k)
            dp1 (k) = delp (i, j, k)
            dp0 (k) = dp1 (k) ! moist air mass * grav

            ! -----------------------------------------------------------------------
            !> - Convert moist mixing ratios to dry mixing ratios.
            ! -----------------------------------------------------------------------

            qvz (k) = qv (i, j, k)
            qlz (k) = ql (i, j, k)
            qrz (k) = qr (i, j, k)
            qgz (k) = qg (i, j, k)

            ! dp1: dry air_mass
            ! dp1 (k) = dp1 (k) * (1. - (qvz (k) + qlz (k) + qrz (k) + qiz (k) + qsz (k) + qgz (k)))
            dp1 (k) = dp1 (k) * (1. - qvz (k)) ! gfs
            omq = dp0 (k) / dp1 (k)

            qvz (k) = qvz (k) * omq
            qlz (k) = qlz (k) * omq
            qrz (k) = qrz (k) * omq
            qiz (k) = qiz (k) * omq
            qsz (k) = qsz (k) * omq
            qgz (k) = qgz (k) * omq

            qa0 (k) = qa (i, j, k)
            qaz (k) = 0.
            dz0 (k) = dz (i, j, k)

            den0 (k) = - dp1 (k) / (grav * dz0 (k)) ! density of dry air
            p1 (k) = den0 (k) * rdgas * t0 (k) ! dry air pressure

            ! -----------------------------------------------------------------------
            ! save a copy of old value for computing tendencies
            ! -----------------------------------------------------------------------

            qv0 (k) = qvz (k)
            ql0 (k) = qlz (k)
            qr0 (k) = qrz (k)
            qi0 (k) = qiz (k)
            qs0 (k) = qsz (k)
            qg0 (k) = qgz (k)

            ! -----------------------------------------------------------------------
            ! for sedi_momentum
            ! -----------------------------------------------------------------------

            m1 (k) = 0.
            u0 (k) = uin (i, j, k)
            v0 (k) = vin (i, j, k)
            u1 (k) = u0 (k)
            v1 (k) = v0 (k)

        enddo

        if (do_sedi_w) then
            do k = ktop, kbot
                w1 (k) = w (i, j, k)
            enddo
        ! Set to -999 if not used
        else
            w1(:) = -999.0
        endif

        ! -----------------------------------------------------------------------
        !> - Calculate cloud condensation nuclei (ccn),
        !! following klein eq. 15
        ! -----------------------------------------------------------------------

        cpaut = c_paut * 0.104 * grav / 1.717e-5

        if (prog_ccn) then
            do k = ktop, kbot
                ! convert # / cc to # / m^3
                ccn (k) = qn (i, j, k) * 1.e6
                c_praut (k) = cpaut * (ccn (k) * rhor) ** (- 1. / 3.)
            enddo
            use_ccn = .false.
        else
            ccn0 = (ccn_l * land (i) + ccn_o * (1. - land (i))) * 1.e6
            if (use_ccn) then
                ! -----------------------------------------------------------------------
                ! ccn is formulted as ccn = ccn_surface * (den / den_surface)
                ! -----------------------------------------------------------------------
                ccn0 = ccn0 * rdgas * tz (kbot) / p1 (kbot)
            endif
            tmp = cpaut * (ccn0 * rhor) ** (- 1. / 3.)
            do k = ktop, kbot
                c_praut (k) = tmp
                ccn (k) = ccn0
            enddo
        endif

        ! -----------------------------------------------------------------------
        !> - Calculate horizontal subgrid variability, which is used in cloud
        !! fraction, relative humidity calculation, evaporation and condensation
        !! processes. Horizontal sub-grid variability is a function of cell area
        !! and land/sea mask:
        !!\n Over land:
        !!\f[
        !! t_{land}=dw_{land}(\frac{A_{r}}{10^{10}})^{0.25}
        !!\f]
        !!\n Over ocean:
        !!\f[
        !! t_{ocean}=dw_{ocean}(\frac{A_{r}}{10^{10}})^{0.25}
        !!\f]
        !! where \f$A_{r}\f$ is cell area. \f$dw_{land}=0.16\f$ and \f$dw_{ocean}=0.10\f$
        !! are base value for sub-grid variability over land and ocean.
        !! The total horizontal sub-grid variability is:
        !!\f[
        !! h_{var}=t_{land}\times fr_{land}+t_{ocean}\times (1-fr_{land})
        !!\f]
        !!\f[
        !! h_{var}=min[0.2,max(0.01,h_{var})]
        !!\f]
        ! total water subgrid deviation in horizontal direction
        ! default area dependent form: use dx ~ 100 km as the base
        ! -----------------------------------------------------------------------

        s_leng = sqrt (sqrt (area1 (i) / 1.e10))
        t_land = dw_land * s_leng
        t_ocean = dw_ocean * s_leng
        h_var = t_land * land (i) + t_ocean * (1. - land (i))
        h_var = min (0.20, max (0.01, h_var))
        ! if (id_var > 0) w_var (i, j) = h_var

        ! -----------------------------------------------------------------------
        !> - Calculate relative humidity increment.
        ! -----------------------------------------------------------------------

        rh_adj = 1. - h_var - rh_inc
        rh_rain = max (0.35, rh_adj - rh_inr) ! rh_inr = 0.25

        ! -----------------------------------------------------------------------
        !> - If requested, call neg_adj() and fix all negative water species.
        ! -----------------------------------------------------------------------

        if (fix_negative) &
            call neg_adj (ktop, kbot, tz, dp1, qvz, qlz, qrz, qiz, qsz, qgz)

        m2_rain (i, :) = 0.
        m2_sol (i, :) = 0.

        !> - Do loop on cloud microphysics sub time step.
        do n = 1, ntimes

            ! -----------------------------------------------------------------------
            !>  - Define air density based on hydrostatical property.
            ! -----------------------------------------------------------------------

            if (p_nonhydro) then
                do k = ktop, kbot
                    dz1 (k) = dz0 (k)
                    den (k) = den0 (k) ! dry air density remains the same
                    denfac (k) = sqrt (sfcrho / den (k))
                enddo
            else
                do k = ktop, kbot
                    dz1 (k) = dz0 (k) * tz (k) / t0 (k) ! hydrostatic balance
                    den (k) = den0 (k) * dz0 (k) / dz1 (k)
                    denfac (k) = sqrt (sfcrho / den (k))
                enddo
            endif

            ! -----------------------------------------------------------------------
            !>  - Call warm_rain() - time-split warm rain processes: 1st pass.
            ! -----------------------------------------------------------------------

            call warm_rain (dt_rain, ktop, kbot, dp1, dz1, tz, qvz, qlz, qrz, qiz, qsz, &
                qgz, den, denfac, ccn, c_praut, rh_rain, vtrz, r1, m1_rain, w1, h_var)

            rain (i) = rain (i) + r1

            do k = ktop, kbot
                m2_rain (i, k) = m2_rain (i, k) + m1_rain (k)
                m1 (k) = m1 (k) + m1_rain (k)
            enddo

            ! -----------------------------------------------------------------------
            !>  - Sedimentation of cloud ice, snow, and graupel.
            ! -----------------------------------------------------------------------

            !>   - Call fall_speed() to calculate the fall velocity of cloud ice, snow
            !!     and graupel.
            call fall_speed (ktop, kbot, den, qsz, qiz, qgz, qlz, tz, vtsz, vtiz, vtgz)

            !>   - Call terminal_fall() to calculate the terminal fall speed.
            call terminal_fall (dts, ktop, kbot, tz, qvz, qlz, qrz, qgz, qsz, qiz, &
                dz1, dp1, den, vtgz, vtsz, vtiz, r1, g1, s1, i1, m1_sol, w1)

            rain (i) = rain (i) + r1 ! from melted snow & ice that reached the ground
            snow (i) = snow (i) + s1
            graupel (i) = graupel (i) + g1
            ice (i) = ice (i) + i1

            ! -----------------------------------------------------------------------
            !>  - Call sedi_heat() to calculate heat transportation during sedimentation.
            ! -----------------------------------------------------------------------

            if (do_sedi_heat) &
                call sedi_heat (ktop, kbot, dp1, m1_sol, dz1, tz, qvz, qlz, qrz, qiz, &
                    qsz, qgz, c_ice)

            ! -----------------------------------------------------------------------
            !>  - Call warm_rain() to - time-split warm rain processes: 2nd pass
            ! -----------------------------------------------------------------------

            call warm_rain (dt_rain, ktop, kbot, dp1, dz1, tz, qvz, qlz, qrz, qiz, qsz, &
                qgz, den, denfac, ccn, c_praut, rh_rain, vtrz, r1, m1_rain, w1, h_var)

            rain (i) = rain (i) + r1

            do k = ktop, kbot
                m2_rain (i, k) = m2_rain (i, k) + m1_rain (k)
                m2_sol (i, k) = m2_sol (i, k) + m1_sol (k)
                m1 (k) = m1 (k) + m1_rain (k) + m1_sol (k)
            enddo

            ! -----------------------------------------------------------------------
            !>  - Call icloud(): ice-phase microphysics
            ! -----------------------------------------------------------------------

            call icloud (ktop, kbot, tz, p1, qvz, qlz, qrz, qiz, qsz, qgz, dp1, den, &
                denfac, vtsz, vtgz, vtrz, qaz, rh_adj, rh_rain, dts, h_var)

        enddo

        ! convert units from Pa*kg/kg to kg/m^2/s
        m2_rain (i, :) = m2_rain (i, :) * rdt * rgrav
        m2_sol (i, :) = m2_sol (i, :) * rdt * rgrav

        ! -----------------------------------------------------------------------
        !> - Calculate momentum transportation during sedimentation.
        ! note: dp1 is dry mass; dp0 is the old moist (total) mass
        ! -----------------------------------------------------------------------

        if (sedi_transport) then
            do k = ktop + 1, kbot
                u1 (k) = (dp0 (k) * u1 (k) + m1 (k - 1) * u1 (k - 1)) / (dp0 (k) + m1 (k - 1))
                v1 (k) = (dp0 (k) * v1 (k) + m1 (k - 1) * v1 (k - 1)) / (dp0 (k) + m1 (k - 1))
                u_dt (i, j, k) = u_dt (i, j, k) + (u1 (k) - u0 (k)) * rdt
                v_dt (i, j, k) = v_dt (i, j, k) + (v1 (k) - v0 (k)) * rdt
            enddo
        endif

        if (do_sedi_w) then
            do k = ktop, kbot
                w (i, j, k) = w1 (k)
            enddo
        endif

        ! -----------------------------------------------------------------------
        !> - Update moist air mass (actually hydrostatic pressure).
        ! convert to dry mixing ratios
        ! -----------------------------------------------------------------------

        do k = ktop, kbot
            omq = dp1 (k) / dp0 (k)
            qv_dt (i, j, k) = qv_dt (i, j, k) + rdt * (qvz (k) - qv0 (k)) * omq
            ql_dt (i, j, k) = ql_dt (i, j, k) + rdt * (qlz (k) - ql0 (k)) * omq
            qr_dt (i, j, k) = qr_dt (i, j, k) + rdt * (qrz (k) - qr0 (k)) * omq
            qi_dt (i, j, k) = qi_dt (i, j, k) + rdt * (qiz (k) - qi0 (k)) * omq
            qs_dt (i, j, k) = qs_dt (i, j, k) + rdt * (qsz (k) - qs0 (k)) * omq
            qg_dt (i, j, k) = qg_dt (i, j, k) + rdt * (qgz (k) - qg0 (k)) * omq
            cvm = c_air + qvz (k) * c_vap + (qrz (k) + qlz (k)) * c_liq + (qiz (k) + qsz (k) + qgz (k)) * c_ice
            pt_dt (i, j, k) = pt_dt (i, j, k) + rdt * (tz (k) - t0 (k)) * cvm / cp_air
        enddo

        ! -----------------------------------------------------------------------
        !> - Update cloud fraction tendency.
        ! -----------------------------------------------------------------------

        do k = ktop, kbot
            if (do_qa) then
                qa_dt (i, j, k) = 0.
            else
                qa_dt (i, j, k) = qa_dt (i, j, k) + rdt * (qaz (k) / real (ntimes) - qa0 (k))
            endif
        enddo

        ! -----------------------------------------------------------------------
        ! fms diagnostics:
        ! -----------------------------------------------------------------------

        ! if (id_cond > 0) then
        ! do k = ktop, kbot ! total condensate
        ! cond (i) = cond (i) + dp1 (k) * (qlz (k) + qrz (k) + qsz (k) + qiz (k) + qgz (k))
        ! enddo
        ! endif
        !
        ! if (id_vtr > 0) then
        ! do k = ktop, kbot
        ! vt_r (i, j, k) = vtrz (k)
        ! enddo
        ! endif
        !
        ! if (id_vts > 0) then
        ! do k = ktop, kbot
        ! vt_s (i, j, k) = vtsz (k)
        ! enddo
        ! endif
        !
        ! if (id_vtg > 0) then
        ! do k = ktop, kbot
        ! vt_g (i, j, k) = vtgz (k)
        ! enddo
        ! endif
        !
        ! if (id_vts > 0) then
        ! do k = ktop, kbot
        ! vt_i (i, j, k) = vtiz (k)
        ! enddo
        ! endif
        !
        ! if (id_droplets > 0) then
        ! do k = ktop, kbot
        ! qn2 (i, j, k) = ccn (k)
        ! enddo
        ! endif

    enddo

end subroutine mpdrv

! -----------------------------------------------------------------------
!>\ingroup mod_gfdl_cloud_mp
!>\brief This subroutine calculates sedimentation of heat.
subroutine sedi_heat (ktop, kbot, dm, m1, dz, tz, qv, ql, qr, qi, qs, qg, cw)

    implicit none

    ! input q fields are dry mixing ratios, and dm is dry air mass

    integer, intent (in) :: ktop, kbot

    real, intent (in), dimension (ktop:kbot) :: dm, m1, dz, qv, ql, qr, qi, qs, qg

    real, intent (inout), dimension (ktop:kbot) :: tz

    real, intent (in) :: cw ! heat capacity

    real, dimension (ktop:kbot) :: dgz, cvn

    real :: tmp

    integer :: k

    do k = ktop, kbot
        dgz (k) = - 0.5 * grav * dz (k) ! > 0
        cvn (k) = dm (k) * (cv_air + qv (k) * cv_vap + (qr (k) + ql (k)) * &
            c_liq + (qi (k) + qs (k) + qg (k)) * c_ice)
    enddo

    ! -----------------------------------------------------------------------
    ! sjl, july 2014
    ! assumption: the ke in the falling condensates is negligible compared to the potential energy
    ! that was unaccounted for. local thermal equilibrium is assumed, and the loss in pe is transformed
    ! into internal energy (to heat the whole grid box)
    ! backward time - implicit upwind transport scheme:
    ! dm here is dry air mass
    ! -----------------------------------------------------------------------

    k = ktop
    tmp = cvn (k) + m1 (k) * cw
    tz (k) = (tmp * tz (k) + m1 (k) * dgz (k)) / tmp

    ! -----------------------------------------------------------------------
    ! implicit algorithm: can't be vectorized
    ! needs an inner i - loop for vectorization
    ! -----------------------------------------------------------------------

    do k = ktop + 1, kbot
        tz (k) = ((cvn (k) + cw * (m1 (k) - m1 (k - 1))) * tz (k) + m1 (k - 1) * &
            cw * tz (k - 1) + dgz (k) * (m1 (k - 1) + m1 (k))) / (cvn (k) + cw * m1 (k))
    enddo

end subroutine sedi_heat

! -----------------------------------------------------------------------
!>\ingroup mod_gfdl_cloud_mp
!> This subroutine includes warm rain cloud microphysics.
!>\section warm_gen GFDL Cloud warm_rain General Algorithm
subroutine warm_rain (dt, ktop, kbot, dp, dz, tz, qv, ql, qr, qi, qs, qg, &
        den, denfac, ccn, c_praut, rh_rain, vtr, r1, m1_rain, w1, h_var)

    implicit none

    integer, intent (in) :: ktop, kbot

    real, intent (in) :: dt ! time step (s)
    real, intent (in) :: rh_rain, h_var

    real, intent (in), dimension (ktop:kbot) :: dp, dz, den
    real, intent (in), dimension (ktop:kbot) :: denfac, ccn, c_praut

    real, intent (inout), dimension (ktop:kbot) :: tz, vtr
    real, intent (inout), dimension (ktop:kbot) :: qv, ql, qr, qi, qs, qg
    real, intent (inout), dimension (ktop:kbot) :: m1_rain, w1

    real, intent (out) :: r1

    real, parameter :: so3 = 7. / 3.

    real, dimension (ktop:kbot) :: dl, dm
    real, dimension (ktop:kbot + 1) :: ze, zt

    real :: sink, dq, qc0, qc
    real :: qden
    real :: zs = 0.
    real :: dt5

    integer :: k

    ! fall velocity constants:

    real, parameter :: vconr = 2503.23638966667
    real, parameter :: normr = 25132741228.7183
    real, parameter :: thr = 1.e-8

    logical :: no_fall

    dt5 = 0.5 * dt

    ! -----------------------------------------------------------------------
    !> - Call check_column() to check if the water species is large enough to fall.
    ! -----------------------------------------------------------------------

    m1_rain (:) = 0.

    call check_column (ktop, kbot, qr, no_fall)

    if (no_fall) then
        vtr (:) = vf_min
        r1 = 0.
    else

        ! -----------------------------------------------------------------------
        !> - Calculate fall speed of rain.
        ! -----------------------------------------------------------------------

        if (const_vr) then
            vtr (:) = vr_fac ! ifs_2016: 4.0
        else
            do k = ktop, kbot
                qden = qr (k) * den (k)
                if (qr (k) < thr) then
                    vtr (k) = vr_min
                else
                    vtr (k) = vr_fac * vconr * sqrt (min (10., sfcrho / den (k))) * &
                        exp (0.2 * log (qden / normr))
                    vtr (k) = min (vr_max, max (vr_min, vtr (k)))
                endif
            enddo
        endif

        ze (kbot + 1) = zs
        do k = kbot, ktop, - 1
            ze (k) = ze (k + 1) - dz (k) ! dz < 0
        enddo

        ! -----------------------------------------------------------------------
        !> - Call revap_racc(), to calculate evaporation and accretion
        !! of rain for the first 1/2 time step.
        ! -----------------------------------------------------------------------

        ! if (.not. fast_sat_adj) &
        call revap_racc (ktop, kbot, dt5, tz, qv, ql, qr, qi, qs, qg, den, denfac, rh_rain, h_var)

        if (do_sedi_w) then
            do k = ktop, kbot
                dm (k) = dp (k) * (1. + qv (k) + ql (k) + qr (k) + qi (k) + qs (k) + qg (k))
            enddo
        endif

        ! -----------------------------------------------------------------------
        !> - Calculate mass flux induced by falling rain
        !! ( use_ppm =.false, call implicit_fall(): time-implicit monotonic fall scheme.)
        ! -----------------------------------------------------------------------

        if (use_ppm) then
            zt (ktop) = ze (ktop)
            do k = ktop + 1, kbot
                zt (k) = ze (k) - dt5 * (vtr (k - 1) + vtr (k))
            enddo
            zt (kbot + 1) = zs - dt * vtr (kbot)

            do k = ktop, kbot
                if (zt (k + 1) >= zt (k)) zt (k + 1) = zt (k) - dz_min
            enddo
            call lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, qr, r1, m1_rain, mono_prof)
        else
            call implicit_fall (dt, ktop, kbot, ze, vtr, dp, qr, r1, m1_rain)
        endif

        ! -----------------------------------------------------------------------
        !  - Calculate vertical velocity transportation during sedimentation.
        ! (do_sedi_w =.true. to turn on vertical motion tranport during sedimentation
        ! .false. by default)
        ! -----------------------------------------------------------------------

        if (do_sedi_w) then
            w1 (ktop) = (dm (ktop) * w1 (ktop) + m1_rain (ktop) * vtr (ktop)) / (dm (ktop) - m1_rain (ktop))
            do k = ktop + 1, kbot
                w1 (k) = (dm (k) * w1 (k) - m1_rain (k - 1) * vtr (k - 1) + m1_rain (k) * vtr (k)) &
                     / (dm (k) + m1_rain (k - 1) - m1_rain (k))
            enddo
        endif

        ! -----------------------------------------------------------------------
        !> - Call sedi_heat() to calculate heat transportation during sedimentation.
        ! -----------------------------------------------------------------------

        if (do_sedi_heat) &
            call sedi_heat (ktop, kbot, dp, m1_rain, dz, tz, qv, ql, qr, qi, qs, qg, c_liq)

        ! -----------------------------------------------------------------------
        !> - Call revap_racc() to calculate evaporation and accretion
        !! of rain for the remaing 1/2 time step.
        ! -----------------------------------------------------------------------

        call revap_racc (ktop, kbot, dt5, tz, qv, ql, qr, qi, qs, qg, den, denfac, rh_rain, h_var)

    endif

    ! -----------------------------------------------------------------------
    !> - Auto-conversion
    !! (assuming linear subgrid vertical distribution of cloud water
    !! following Lin et al. (1994) \cite lin_et_al_1994.)
    ! -----------------------------------------------------------------------

    if (irain_f /= 0) then

        ! -----------------------------------------------------------------------
        ! no subgrid varaibility
        ! -----------------------------------------------------------------------

        do k = ktop, kbot
            qc0 = fac_rc * ccn (k)
            if (tz (k) > t_wfr) then
                if (use_ccn) then
                    ! -----------------------------------------------------------------------
                    ! ccn is formulted as ccn = ccn_surface * (den / den_surface)
                    ! -----------------------------------------------------------------------
                    qc = qc0
                else
                    qc = qc0 / den (k)
                endif
                dq = ql (k) - qc
                if (dq > 0.) then
                    sink = min (dq, dt * c_praut (k) * den (k) * exp (so3 * log (ql (k))))
                    ql (k) = ql (k) - sink
                    qr (k) = qr (k) + sink
                endif
            endif
        enddo

    else

        ! -----------------------------------------------------------------------
        !>  - Call linear_prof() to calculate vertical subgrid variability of cloud water.
        ! -----------------------------------------------------------------------

        call linear_prof (kbot - ktop + 1, ql (ktop), dl (ktop), z_slope_liq, h_var)

        do k = ktop, kbot
            qc0 = fac_rc * ccn (k)
            if (tz (k) > t_wfr + dt_fr) then
                dl (k) = min (max (1.e-6, dl (k)), 0.5 * ql (k))
                ! --------------------------------------------------------------------
                ! as in klein's gfdl am2 stratiform scheme (with subgrid variations)
                ! --------------------------------------------------------------------
                if (use_ccn) then
                    ! --------------------------------------------------------------------
                    ! ccn is formulted as ccn = ccn_surface * (den / den_surface)
                    ! --------------------------------------------------------------------
                    qc = qc0
                else
                    qc = qc0 / den (k)
                endif
                dq = 0.5 * (ql (k) + dl (k) - qc)
                ! --------------------------------------------------------------------
                ! dq = dl if qc == q_minus = ql - dl
                ! dq = 0 if qc == q_plus = ql + dl
                ! --------------------------------------------------------------------
                if (dq > 0.) then ! q_plus > qc
                    ! --------------------------------------------------------------------
                    !>  - Revised continuous form: linearly decays (with subgrid dl) to zero at qc == ql + dl
                    ! --------------------------------------------------------------------
                    sink = min (1., dq / dl (k)) * dt * c_praut (k) * den (k) * exp (so3 * log (ql (k)))
                    ql (k) = ql (k) - sink
                    qr (k) = qr (k) + sink
                endif
            endif
        enddo
    endif

end subroutine warm_rain

! -----------------------------------------------------------------------
!>\ingroup mod_gfdl_cloud_mp
!> This subroutine calculates evaporation of rain and accretion of rain.
!!\section gen_ravap GFDL Cloud revap_racc General Algorithm
subroutine revap_racc (ktop, kbot, dt, tz, qv, ql, qr, qi, qs, qg, den, denfac, rh_rain, h_var)

    implicit none

    integer, intent (in) :: ktop, kbot

    real, intent (in) :: dt ! time step (s)
    real, intent (in) :: rh_rain, h_var

    real, intent (in), dimension (ktop:kbot) :: den, denfac

    real, intent (inout), dimension (ktop:kbot) :: tz, qv, qr, ql, qi, qs, qg

    real, dimension (ktop:kbot) :: lhl, cvm, q_liq, q_sol, lcpk

    real :: dqv, qsat, dqsdt, evap, t2, qden, q_plus, q_minus, sink
    real :: qpz, dq, dqh, tin

    integer :: k

    do k = ktop, kbot

        if (tz (k) > t_wfr .and. qr (k) > qrmin) then

            ! -----------------------------------------------------------------------
            !> - Define heat capacity and latent heat coefficient.
            ! -----------------------------------------------------------------------

            lhl (k) = lv00 + d0_vap * tz (k)
            q_liq (k) = ql (k) + qr (k)
            q_sol (k) = qi (k) + qs (k) + qg (k)
            cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
            lcpk (k) = lhl (k) / cvm (k)

            tin = tz (k) - lcpk (k) * ql (k) ! presence of clouds suppresses the rain evap
            qpz = qv (k) + ql (k)
            qsat = wqs2 (tin, den (k), dqsdt)
            dqh = max (ql (k), h_var * max (qpz, qcmin))
            dqh = min (dqh, 0.2 * qpz) ! new limiter
            dqv = qsat - qv (k) ! use this to prevent super - sat the gird box
            q_minus = qpz - dqh
            q_plus = qpz + dqh

            ! -----------------------------------------------------------------------
            !> - qsat must be > q_minus to activate evaporation;
            !! qsat must be < q_plus to activate accretion
            ! -----------------------------------------------------------------------

            ! -----------------------------------------------------------------------
            !> - rain evaporation (\f$evap\f$)
            ! -----------------------------------------------------------------------

            if (dqv > qvmin .and. qsat > q_minus) then
                if (qsat > q_plus) then
                    dq = qsat - qpz
                else
                    ! -----------------------------------------------------------------------
                    ! q_minus < qsat < q_plus
                    ! dq == dqh if qsat == q_minus
                    ! -----------------------------------------------------------------------
                    dq = 0.25 * (q_minus - qsat) ** 2 / dqh
                endif
                qden = qr (k) * den (k)
                t2 = tin * tin
                evap = crevp (1) * t2 * dq * (crevp (2) * sqrt (qden) + crevp (3) * &
                    exp (0.725 * log (qden))) / (crevp (4) * t2 + crevp (5) * qsat * den (k))
                evap = min (qr (k), dt * evap, dqv / (1. + lcpk (k) * dqsdt))
                ! -----------------------------------------------------------------------
                ! alternative minimum evap in dry environmental air
                ! sink = min (qr (k), dim (rh_rain * qsat, qv (k)) / (1. + lcpk (k) * dqsdt))
                ! evap = max (evap, sink)
                ! -----------------------------------------------------------------------
                qr (k) = qr (k) - evap
                qv (k) = qv (k) + evap
                q_liq (k) = q_liq (k) - evap
                cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
                tz (k) = tz (k) - evap * lhl (k) / cvm (k)
            endif

            ! -----------------------------------------------------------------------
            !> - accretion: \f$P_{racc}\f$
            ! -----------------------------------------------------------------------

            ! if (qr (k) > qrmin .and. ql (k) > 1.e-7 .and. qsat < q_plus) then
            if (qr (k) > qrmin .and. ql (k) > 1.e-6 .and. qsat < q_minus) then
                sink = dt * denfac (k) * cracw * exp (0.95 * log (qr (k) * den (k)))
                sink = sink / (1. + sink) * ql (k)
                ql (k) = ql (k) - sink
                qr (k) = qr (k) + sink
            endif

        endif ! warm - rain
    enddo

end subroutine revap_racc

! -----------------------------------------------------------------------
!>\ingroup mod_gfdl_cloud_mp
!> Definition of vertical subgrid variability
!! used for cloud ice and cloud water autoconversion.
! qi -- > ql & ql -- > qr
! edges: qe == qbar + / - dm
!>\section gen_linear GFDL cloud linear_prof General Algorithm
subroutine linear_prof (km, q, dm, z_var, h_var)

    implicit none

    integer, intent (in) :: km

    real, intent (in) :: q (km), h_var

    real, intent (out) :: dm (km)

    logical, intent (in) :: z_var

    real :: dq (km)

    integer :: k

    if (z_var) then
        do k = 2, km
            dq (k) = 0.5 * (q (k) - q (k - 1))
        enddo
        dm (1) = 0.

        ! -----------------------------------------------------------------------
        !> - Use twice the strength of the positive definiteness limiter (Lin et al.(1994)
        !! \cite lin_et_al_1994).
        ! -----------------------------------------------------------------------

        do k = 2, km - 1
            dm (k) = 0.5 * min (abs (dq (k) + dq (k + 1)), 0.5 * q (k))
            if (dq (k) * dq (k + 1) <= 0.) then
                if (dq (k) > 0.) then ! local max
                    dm (k) = min (dm (k), dq (k), - dq (k + 1))
                else
                    dm (k) = 0.
                endif
            endif
        enddo
        dm (km) = 0.

        ! -----------------------------------------------------------------------
        !> - Impose the background horizontal variability (\f$h_{var}\f$) that
        !! is proportional to the value itself.
        ! -----------------------------------------------------------------------

        do k = 1, km
            dm (k) = max (dm (k), qvmin, h_var * q (k))
        enddo
    else
        do k = 1, km
            dm (k) = max (qvmin, h_var * q (k))
        enddo
    endif

end subroutine linear_prof

! =======================================================================
!>\ingroup mod_gfdl_cloud_mp
!> This subroutine includes cloud ice microphysics processes.
!>\author Shian-Jiann Lin, GFDL
!!
!! This scheme is featured with:
!! - bulk cloud microphysics
!! - processes splitting with some un-split sub-grouping
!! - time implicit (when possible) accretion and autoconversion
!>\section det_icloud GFDL icloud Detailed Algorithm
subroutine icloud (ktop, kbot, tzk, p1, qvk, qlk, qrk, qik, qsk, qgk, dp1, &
        den, denfac, vts, vtg, vtr, qak, rh_adj, rh_rain, dts, h_var)

    implicit none

    integer, intent (in) :: ktop, kbot

    real, intent (in), dimension (ktop:kbot) :: p1, dp1, den, denfac, vts, vtg, vtr

    real, intent (inout), dimension (ktop:kbot) :: tzk, qvk, qlk, qrk, qik, qsk, qgk, qak

    real, intent (in) :: rh_adj, rh_rain, dts, h_var

    real, dimension (ktop:kbot) :: lcpk, icpk, tcpk, di, lhl, lhi
    real, dimension (ktop:kbot) :: cvm, q_liq, q_sol

    real :: rdts, fac_g2v, fac_v2g, fac_i2s, fac_imlt
    real :: tz, qv, ql, qr, qi, qs, qg, melt
    real :: pracs, psacw, pgacw, psacr, pgacr, pgaci, praci, psaci
    real :: pgmlt, psmlt, pgfr, pgaut, psaut, pgsub
    real :: tc, tsq, dqs0, qden, qim, qsm
    real :: dt5, factor, sink, qi_crt
    real :: tmp, qsw, qsi, dqsdt, dq
    real :: dtmp, qc, q_plus, q_minus

    integer :: k

    dt5 = 0.5 * dts

    rdts = 1. / dts

    ! -----------------------------------------------------------------------
    !> - Define conversion scalar/factor.
    ! -----------------------------------------------------------------------

    fac_i2s = 1. - exp (- dts / tau_i2s)
    fac_g2v = 1. - exp (- dts / tau_g2v)
    fac_v2g = 1. - exp (- dts / tau_v2g)

    fac_imlt = 1. - exp (- dt5 / tau_imlt)

    ! -----------------------------------------------------------------------
    !> - Define heat capacity and latend heat coefficient.
    ! -----------------------------------------------------------------------

    do k = ktop, kbot
        lhi (k) = li00 + dc_ice * tzk (k)
        q_liq (k) = qlk (k) + qrk (k)
        q_sol (k) = qik (k) + qsk (k) + qgk (k)
        cvm (k) = c_air + qvk (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
        icpk (k) = lhi (k) / cvm (k)
    enddo

    ! -----------------------------------------------------------------------
    ! sources of cloud ice: pihom, cold rain, and the sat_adj
    ! (initiation plus deposition)
    ! sources of snow: cold rain, auto conversion + accretion (from cloud ice)
    ! sat_adj (deposition; requires pre - existing snow) ; initial snow comes from auto conversion
    ! -----------------------------------------------------------------------

    do k = ktop, kbot
        if (tzk (k) > tice .and. qik (k) > qcmin) then

            ! -----------------------------------------------------------------------
            !> - Calculate \f$P_{imlt}\f$: instant melting of cloud ice.
            ! -----------------------------------------------------------------------

            melt = min (qik (k), fac_imlt * (tzk (k) - tice) / icpk (k))
            tmp = min (melt, dim (ql_mlt, qlk (k))) ! max ql amount
            qlk (k) = qlk (k) + tmp
            qrk (k) = qrk (k) + melt - tmp
            qik (k) = qik (k) - melt
            q_liq (k) = q_liq (k) + melt
            q_sol (k) = q_sol (k) - melt
            cvm (k) = c_air + qvk (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
            tzk (k) = tzk (k) - melt * lhi (k) / cvm (k)

        elseif (tzk (k) < t_wfr .and. qlk (k) > qcmin) then

            ! -----------------------------------------------------------------------
            !> - Calculate \f$P_{ihom}\f$: homogeneous freezing of cloud water into cloud ice.
            !! This is the 1st occurance of liquid water freezing in the split MP process.
            ! -----------------------------------------------------------------------

            dtmp = t_wfr - tzk (k)
            factor = min (1., dtmp / dt_fr)
            sink = min (qlk (k) * factor, dtmp / icpk (k))
            qi_crt = qi_gen * min (qi_lim, 0.1 * (tice - tzk (k))) / den (k)
            tmp = min (sink, dim (qi_crt, qik (k)))
            qlk (k) = qlk (k) - sink
            qsk (k) = qsk (k) + sink - tmp
            qik (k) = qik (k) + tmp
            q_liq (k) = q_liq (k) - sink
            q_sol (k) = q_sol (k) + sink
            cvm (k) = c_air + qvk (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
            tzk (k) = tzk (k) + sink * lhi (k) / cvm (k)

        endif
    enddo

    ! -----------------------------------------------------------------------
    !> - Call linear_prof() to calculate vertical subgrid variability of cloud ice.
    ! -----------------------------------------------------------------------

    call linear_prof (kbot - ktop + 1, qik (ktop), di (ktop), z_slope_ice, h_var)

    ! -----------------------------------------------------------------------
    !> - Update capacity heat and latend heat coefficient.
    ! -----------------------------------------------------------------------

    do k = ktop, kbot
        lhl (k) = lv00 + d0_vap * tzk (k)
        lhi (k) = li00 + dc_ice * tzk (k)
        lcpk (k) = lhl (k) / cvm (k)
        icpk (k) = lhi (k) / cvm (k)
        tcpk (k) = lcpk (k) + icpk (k)
    enddo

    do k = ktop, kbot

        ! -----------------------------------------------------------------------
        ! do nothing above p_min
        ! -----------------------------------------------------------------------

        if (p1 (k) < p_min) cycle

        tz = tzk (k)
        qv = qvk (k)
        ql = qlk (k)
        qi = qik (k)
        qr = qrk (k)
        qs = qsk (k)
        qg = qgk (k)

        pgacr = 0.
        pgacw = 0.
        tc = tz - tice

        if (tc .ge. 0.) then

            ! -----------------------------------------------------------------------
            !> - Melting of snow:
            ! -----------------------------------------------------------------------

            dqs0 = ces0 / p1 (k) - qv

            if (qs > qcmin) then

                ! -----------------------------------------------------------------------
                !>  - \f$P_{sacw}\f$: accretion of cloud water by snow
                ! only rate is used (for snow melt) since tc > 0.
                ! -----------------------------------------------------------------------

                if (ql > qrmin) then
                    factor = denfac (k) * csacw * exp (0.8125 * log (qs * den (k)))
                    psacw = factor / (1. + dts * factor) * ql ! rate
                else
                    psacw = 0.
                endif

                ! -----------------------------------------------------------------------
                !>  - \f$P_{sacr}\f$: accretion of rain by melted snow
                !>  - \f$P_{racs}\f$: accretion of snow by rain
                ! -----------------------------------------------------------------------

                if (qr > qrmin) then
                    psacr = min (acr3d (vts (k), vtr (k), qr, qs, csacr, acco (1, 2), &
                        den (k)), qr * rdts)
                    pracs = acr3d (vtr (k), vts (k), qs, qr, cracs, acco (1, 1), den (k))
                else
                    psacr = 0.
                    pracs = 0.
                endif

                ! -----------------------------------------------------------------------
                !>  - Total snow sink:
                !! \f$P_{smlt}\f$: snow melt (smlt(); due to rain accretion)
                ! -----------------------------------------------------------------------

                psmlt = max (0., smlt (tc, dqs0, qs * den (k), psacw, psacr, csmlt, &
                    den (k), denfac (k)))
                sink = min (qs, dts * (psmlt + pracs), tc / icpk (k))
                qs = qs - sink
                ! sjl, 20170321:
                tmp = min (sink, dim (qs_mlt, ql)) ! max ql due to snow melt
                ql = ql + tmp
                qr = qr + sink - tmp
                ! qr = qr + sink
                ! sjl, 20170321:
                q_liq (k) = q_liq (k) + sink
                q_sol (k) = q_sol (k) - sink
                cvm (k) = c_air + qv * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
                tz = tz - sink * lhi (k) / cvm (k)
                tc = tz - tice

            endif

            ! -----------------------------------------------------------------------
            !> - Update capacity heat and latend heat coefficient.
            ! -----------------------------------------------------------------------

            lhi (k) = li00 + dc_ice * tz
            icpk (k) = lhi (k) / cvm (k)

            ! -----------------------------------------------------------------------
            !> - Melting of graupel:
            ! -----------------------------------------------------------------------

            if (qg > qcmin .and. tc > 0.) then

                ! -----------------------------------------------------------------------
                !>  - \f$P_{gacr}\f$: accretion of rain by graupel
                ! -----------------------------------------------------------------------

                if (qr > qrmin) &
                    pgacr = min (acr3d (vtg (k), vtr (k), qr, qg, cgacr, acco (1, 3), &
                        den (k)), rdts * qr)

                ! -----------------------------------------------------------------------
                !>  - \f$P_{gacw}\f$: accretion of cloud water by graupel
                ! -----------------------------------------------------------------------

                qden = qg * den (k)
                if (ql > qrmin) then
                    factor = cgacw * qden / sqrt (den (k) * sqrt (sqrt (qden)))
                    pgacw = factor / (1. + dts * factor) * ql ! rate
                endif

                ! -----------------------------------------------------------------------
                !>  - \f$P_{gmlt}\f$: graupel melt (gmlt())
                ! -----------------------------------------------------------------------

                pgmlt = dts * gmlt (tc, dqs0, qden, pgacw, pgacr, cgmlt, den (k))
                pgmlt = min (max (0., pgmlt), qg, tc / icpk (k))
                qg = qg - pgmlt
                qr = qr + pgmlt
                q_liq (k) = q_liq (k) + pgmlt
                q_sol (k) = q_sol (k) - pgmlt
                cvm (k) = c_air + qv * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
                tz = tz - pgmlt * lhi (k) / cvm (k)

            endif

        else

            ! -----------------------------------------------------------------------
            !> - Cloud ice processes:
            ! -----------------------------------------------------------------------

            ! -----------------------------------------------------------------------
            !>  - \f$P_{saci}\f$: accretion of cloud ice by snow
            ! -----------------------------------------------------------------------

            if (qi > 3.e-7) then ! cloud ice sink terms

                if (qs > 1.e-7) then
                    ! -----------------------------------------------------------------------
                    ! sjl added (following lin eq. 23) the temperature dependency
                    ! to reduce accretion, use esi = exp (0.05 * tc) as in hong et al 2004
                    ! -----------------------------------------------------------------------
                    factor = dts * denfac (k) * csaci * exp (0.05 * tc + 0.8125 * log (qs * den (k)))
                    psaci = factor / (1. + factor) * qi
                else
                    psaci = 0.
                endif

                ! -----------------------------------------------------------------------
                !>  - \f$P_{saut}\f$: autoconversion: cloud ice \f$\rightarrow\f$ snow.
                !!\n similar to Lin et al.(1983) \cite lin_et_al_1983 eq. 21 solved implicitly;
                !! threshold from wsm6 scheme, Hong et al. (2004) \cite hong_et_al_2004,
                !! eq (13) : qi0_crt ~0.8e-4.
                ! -----------------------------------------------------------------------
                if (qi0_crt < 0.) then
                    qim = - qi0_crt
                else
                    qim = qi0_crt / den (k)
                endif
                ! -----------------------------------------------------------------------
                ! assuming linear subgrid vertical distribution of cloud ice
                ! the mismatch computation following lin et al. 1994, mwr
                ! -----------------------------------------------------------------------

                if (const_vi) then
                    tmp = fac_i2s
                else
                    tmp = fac_i2s * exp (0.025 * tc)
                endif

                di (k) = max (di (k), qrmin)
                q_plus = qi + di (k)
                if (q_plus > (qim + qrmin)) then
                    if (qim > (qi - di (k))) then
                        dq = (0.25 * (q_plus - qim) ** 2) / di (k)
                    else
                        dq = qi - qim
                    endif
                    psaut = tmp * dq
                else
                    psaut = 0.
                endif
                ! -----------------------------------------------------------------------
                ! sink is no greater than 75% of qi
                ! -----------------------------------------------------------------------
                sink = min (0.75 * qi, psaci + psaut)
                qi = qi - sink
                qs = qs + sink

                ! -----------------------------------------------------------------------
                !>  - \f$P_{gaci}\f$: accretion of cloud ice by graupel
                ! -----------------------------------------------------------------------

                if (qg > 1.e-6) then
                    ! -----------------------------------------------------------------------
                    ! factor = dts * cgaci / sqrt (den (k)) * exp (0.05 * tc + 0.875 * log (qg * den (k)))
                    ! simplified form: remove temp dependency & set the exponent "0.875" -- > 1
                    ! -----------------------------------------------------------------------
                    factor = dts * cgaci * sqrt (den (k)) * qg
                    pgaci = factor / (1. + factor) * qi
                    qi = qi - pgaci
                    qg = qg + pgaci
                endif

            endif

            ! -----------------------------------------------------------------------
            !> - Cold-rain processes:
            ! -----------------------------------------------------------------------

            ! -----------------------------------------------------------------------
            ! rain to ice, snow, graupel processes:
            ! -----------------------------------------------------------------------

            tc = tz - tice

            if (qr > 1.e-7 .and. tc < 0.) then

                ! -----------------------------------------------------------------------
                ! * sink * terms to qr: psacr + pgfr
                ! source terms to qs: psacr
                ! source terms to qg: pgfr
                ! -----------------------------------------------------------------------

                ! -----------------------------------------------------------------------
                !>  - \f$P_{sacr}\f$: accretion of rain by snow (acr3d())
                ! -----------------------------------------------------------------------

                if (qs > 1.e-7) then ! if snow exists
                    psacr = dts * acr3d (vts (k), vtr (k), qr, qs, csacr, acco (1, 2), den (k))
                else
                    psacr = 0.
                endif

                ! -----------------------------------------------------------------------
                !>  - \f$P_{gfr}\f$: rain freezing \f$\rightarrow\f$ graupel
                ! -----------------------------------------------------------------------

                pgfr = dts * cgfr (1) / den (k) * (exp (- cgfr (2) * tc) - 1.) * &
                    exp (1.75 * log (qr * den (k)))

                ! -----------------------------------------------------------------------
                !>  - Calculate total sink to \f$q_r\f$ :
                !!\n  Sink terms to \f$q_r\f$:  \f$P_{sacr}+P_{gfr}\f$
                !!\n  source term to \f$q_s\f$: \f$P_{sacr}\f$
                !!\n  source term to \f$q_g\f$: \f$P_{gfr}\f$
                ! -----------------------------------------------------------------------

                sink = psacr + pgfr
                factor = min (sink, qr, - tc / icpk (k)) / max (sink, qrmin)

                psacr = factor * psacr
                pgfr = factor * pgfr

                sink = psacr + pgfr
                qr = qr - sink
                qs = qs + psacr
                qg = qg + pgfr
                q_liq (k) = q_liq (k) - sink
                q_sol (k) = q_sol (k) + sink
                cvm (k) = c_air + qv * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
                tz = tz + sink * lhi (k) / cvm (k)

            endif

            ! -----------------------------------------------------------------------
            !> - Update capacity heat and latend heat coefficient.
            ! -----------------------------------------------------------------------

            lhi (k) = li00 + dc_ice * tz
            icpk (k) = lhi (k) / cvm (k)

            ! -----------------------------------------------------------------------
            !> - Graupel production terms:
            ! -----------------------------------------------------------------------

            if (qs > 1.e-7) then

                ! -----------------------------------------------------------------------
                !>  - accretion: snow \f$\rightarrow\f$ graupel (acr3d())
                ! -----------------------------------------------------------------------

                if (qg > qrmin) then
                    sink = dts * acr3d (vtg (k), vts (k), qs, qg, cgacs, acco (1, 4), den (k))
                else
                    sink = 0.
                endif

                ! -----------------------------------------------------------------------
                !>  - autoconversion: snow \f$\rightarrow\f$ graupel
                ! -----------------------------------------------------------------------

                qsm = qs0_crt / den (k)
                if (qs > qsm) then
                    factor = dts * 1.e-3 * exp (0.09 * (tz - tice))
                    sink = sink + factor / (1. + factor) * (qs - qsm)
                endif
                sink = min (qs, sink)
                qs = qs - sink
                qg = qg + sink

            endif ! snow existed

            if (qg > 1.e-7 .and. tz < tice0) then

                ! -----------------------------------------------------------------------
                !>  - \f$P_{gacw}\f$: accretion of cloud water by graupel
                ! -----------------------------------------------------------------------

                if (ql > 1.e-6) then
                    qden = qg * den (k)
                    factor = dts * cgacw * qden / sqrt (den (k) * sqrt (sqrt (qden)))
                    pgacw = factor / (1. + factor) * ql
                else
                    pgacw = 0.
                endif

                ! -----------------------------------------------------------------------
                !>  - \f$P_{gacr}\f$: accretion of rain by graupel (acr3d())
                ! -----------------------------------------------------------------------

                if (qr > 1.e-6) then
                    pgacr = min (dts * acr3d (vtg (k), vtr (k), qr, qg, cgacr, acco (1, 3), &
                        den (k)), qr)
                else
                    pgacr = 0.
                endif

                sink = pgacr + pgacw
                factor = min (sink, dim (tice, tz) / icpk (k)) / max (sink, qrmin)
                pgacr = factor * pgacr
                pgacw = factor * pgacw

                sink = pgacr + pgacw
                qg = qg + sink
                qr = qr - pgacr
                ql = ql - pgacw
                q_liq (k) = q_liq (k) - sink
                q_sol (k) = q_sol (k) + sink
                cvm (k) = c_air + qv * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
                tz = tz + sink * lhi (k) / cvm (k)

            endif

        endif

        tzk (k) = tz
        qvk (k) = qv
        qlk (k) = ql
        qik (k) = qi
        qrk (k) = qr
        qsk (k) = qs
        qgk (k) = qg

    enddo

    ! -----------------------------------------------------------------------
    !> - Call subgrid_z_proc() for subgrid cloud microphysics.
    ! -----------------------------------------------------------------------

    call subgrid_z_proc (ktop, kbot, p1, den, denfac, dts, rh_adj, tzk, qvk, &
        qlk, qrk, qik, qsk, qgk, qak, h_var, rh_rain)

end subroutine icloud

! =======================================================================
!>\ingroup mod_gfdl_cloud_mp
!> This subroutine calculates temperature sentive high vertical resolution processes.
!>\section gen_subz GFDL Cloud subgrid_z_proc General Algorithm
subroutine subgrid_z_proc (ktop, kbot, p1, den, denfac, dts, rh_adj, tz, qv, &
    ql, qr, qi, qs, qg, qa, h_var, rh_rain)

    implicit none

    integer, intent (in) :: ktop, kbot

    real, intent (in), dimension (ktop:kbot) :: p1, den, denfac

    real, intent (in) :: dts, rh_adj, h_var, rh_rain

    real, intent (inout), dimension (ktop:kbot) :: tz, qv, ql, qr, qi, qs, qg, qa

    real, dimension (ktop:kbot) :: lcpk, icpk, tcpk, tcp3, lhl, lhi
    real, dimension (ktop:kbot) :: cvm, q_liq, q_sol, q_cond

    real :: fac_v2l, fac_l2v

    real :: pidep, qi_crt

    ! -----------------------------------------------------------------------
    ! qstar over water may be accurate only down to - 80 deg c with ~10% uncertainty
    ! must not be too large to allow psc
    ! -----------------------------------------------------------------------

    real :: rh, rqi, tin, qsw, qsi, qpz, qstar
    real :: dqsdt, dwsdt, dq, dq0, factor, tmp
    real :: q_plus, q_minus, dt_evap, dt_pisub
    real :: evap, sink, tc, pisub, q_adj, dtmp
    real :: pssub, pgsub, tsq, qden, fac_g2v, fac_v2g

    integer :: k

    if (fast_sat_adj) then
        dt_evap = 0.5 * dts
    else
        dt_evap = dts
    endif

    ! -----------------------------------------------------------------------
    !> - Define conversion scalar/factor.
    ! -----------------------------------------------------------------------

    fac_v2l = 1. - exp (- dt_evap / tau_v2l)
    fac_l2v = 1. - exp (- dt_evap / tau_l2v)

    fac_g2v = 1. - exp (- dts / tau_g2v)
    fac_v2g = 1. - exp (- dts / tau_v2g)

    ! -----------------------------------------------------------------------
    !> - Define heat capacity and latend heat coefficient.
    ! -----------------------------------------------------------------------

    do k = ktop, kbot
        lhl (k) = lv00 + d0_vap * tz (k)
        lhi (k) = li00 + dc_ice * tz (k)
        q_liq (k) = ql (k) + qr (k)
        q_sol (k) = qi (k) + qs (k) + qg (k)
        cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
        lcpk (k) = lhl (k) / cvm (k)
        icpk (k) = lhi (k) / cvm (k)
        tcpk (k) = lcpk (k) + icpk (k)
        tcp3 (k) = lcpk (k) + icpk (k) * min (1., dim (tice, tz (k)) / (tice - t_wfr))
    enddo

    do k = ktop, kbot

        if (p1 (k) < p_min) cycle

        ! -----------------------------------------------------------------------
        !> - Instant deposit all water vapor to cloud ice when temperature is super low.
        ! -----------------------------------------------------------------------

        if (tz (k) < t_min) then
            sink = dim (qv (k), 1.e-7)
            qv (k) = qv (k) - sink
            qi (k) = qi (k) + sink
            q_sol (k) = q_sol (k) + sink
            cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
            tz (k) = tz (k) + sink * (lhl (k) + lhi (k)) / cvm (k)
            if (.not. do_qa) qa (k) = qa (k) + 1. ! air fully saturated; 100 % cloud cover
            cycle
        endif

        ! -----------------------------------------------------------------------
        !> - Update heat capacity and latend heat coefficient.
        ! -----------------------------------------------------------------------

        lhl (k) = lv00 + d0_vap * tz (k)
        lhi (k) = li00 + dc_ice * tz (k)
        lcpk (k) = lhl (k) / cvm (k)
        icpk (k) = lhi (k) / cvm (k)
        tcpk (k) = lcpk (k) + icpk (k)
        tcp3 (k) = lcpk (k) + icpk (k) * min (1., dim (tice, tz (k)) / (tice - t_wfr))

        ! -----------------------------------------------------------------------
        !> - Instant evaporation/sublimation of all clouds if rh < rh_adj \f$\rightarrow\f$ cloud free.
        ! -----------------------------------------------------------------------

        qpz = qv (k) + ql (k) + qi (k)
        tin = tz (k) - (lhl (k) * (ql (k) + qi (k)) + lhi (k) * qi (k)) / (c_air + &
            qpz * c_vap + qr (k) * c_liq + (qs (k) + qg (k)) * c_ice)
        if (tin > t_sub + 6.) then
            rh = qpz / iqs1 (tin, den (k))
            if (rh < rh_adj) then ! qpz / rh_adj < qs
                tz (k) = tin
                qv (k) = qpz
                ql (k) = 0.
                qi (k) = 0.
                cycle ! cloud free
            endif
        endif

        ! -----------------------------------------------------------------------
        !> - cloud water \f$\Leftrightarrow\f$ vapor adjustment:
        ! -----------------------------------------------------------------------

        qsw = wqs2 (tz (k), den (k), dwsdt)
        dq0 = qsw - qv (k)
        if (dq0 > 0.) then
            ! SJL 20170703 added ql factor to prevent the situation of high ql and low RH
            ! factor = min (1., fac_l2v * sqrt (max (0., ql (k)) / 1.e-5) * 10. * dq0 / qsw)
            ! factor = fac_l2v
            ! factor = 1
            factor = min (1., fac_l2v * (10. * dq0 / qsw)) ! the rh dependent factor = 1 at 90%
            evap = min (ql (k), factor * dq0 / (1. + tcp3 (k) * dwsdt))
        else ! condensate all excess vapor into cloud water
            ! -----------------------------------------------------------------------
            ! evap = fac_v2l * dq0 / (1. + tcp3 (k) * dwsdt)
            ! sjl, 20161108
            ! -----------------------------------------------------------------------
            evap = dq0 / (1. + tcp3 (k) * dwsdt)
        endif
        qv (k) = qv (k) + evap
        ql (k) = ql (k) - evap
        q_liq (k) = q_liq (k) - evap
        cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
        tz (k) = tz (k) - evap * lhl (k) / cvm (k)

        ! -----------------------------------------------------------------------
        !> - Update heat capacity and latend heat coefficient.
        ! -----------------------------------------------------------------------

        lhi (k) = li00 + dc_ice * tz (k)
        icpk (k) = lhi (k) / cvm (k)

        ! -----------------------------------------------------------------------
        !> - Enforce complete freezing below \f$-48^oC\f$.
        ! -----------------------------------------------------------------------

        dtmp = t_wfr - tz (k) ! [ - 40, - 48]
        if (dtmp > 0. .and. ql (k) > qcmin) then
            sink = min (ql (k), ql (k) * dtmp * 0.125, dtmp / icpk (k))
            ql (k) = ql (k) - sink
            qi (k) = qi (k) + sink
            q_liq (k) = q_liq (k) - sink
            q_sol (k) = q_sol (k) + sink
            cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
            tz (k) = tz (k) + sink * lhi (k) / cvm (k)
        endif

        ! -----------------------------------------------------------------------
        !> - Update heat capacity and latend heat coefficient.
        ! -----------------------------------------------------------------------

        lhi (k) = li00 + dc_ice * tz (k)
        icpk (k) = lhi (k) / cvm (k)

        ! -----------------------------------------------------------------------
        !> - Apply Bigg mechanism.
        ! -----------------------------------------------------------------------

        if (fast_sat_adj) then
            dt_pisub = 0.5 * dts
        else
            dt_pisub = dts
            tc = tice - tz (k)
            if (ql (k) > qrmin .and. tc > 0.) then
                sink = 3.3333e-10 * dts * (exp (0.66 * tc) - 1.) * den (k) * ql (k) * ql (k)
                sink = min (ql (k), tc / icpk (k), sink)
                ql (k) = ql (k) - sink
                qi (k) = qi (k) + sink
                q_liq (k) = q_liq (k) - sink
                q_sol (k) = q_sol (k) + sink
                cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
                tz (k) = tz (k) + sink * lhi (k) / cvm (k)
            endif ! significant ql existed
        endif

        ! -----------------------------------------------------------------------
        !> - Update capacity heat and latend heat coefficient.
        ! -----------------------------------------------------------------------

        lhl (k) = lv00 + d0_vap * tz (k)
        lhi (k) = li00 + dc_ice * tz (k)
        lcpk (k) = lhl (k) / cvm (k)
        icpk (k) = lhi (k) / cvm (k)
        tcpk (k) = lcpk (k) + icpk (k)

        ! -----------------------------------------------------------------------
        !> - Sublimation/deposition of ice.
        ! -----------------------------------------------------------------------

        if (tz (k) < tice) then
            qsi = iqs2 (tz (k), den (k), dqsdt)
            dq = qv (k) - qsi
            sink = dq / (1. + tcpk (k) * dqsdt)
            if (qi (k) > qrmin) then
                ! eq 9, hong et al. 2004, mwr
                ! for a and b, see dudhia 1989: page 3103 eq (b7) and (b8)
                pidep = dt_pisub * dq * 349138.78 * exp (0.875 * log (qi (k) * den (k))) &
                     / (qsi * den (k) * lat2 / (0.0243 * rvgas * tz (k) ** 2) + 4.42478e4)
            else
                pidep = 0.
            endif
            if (dq > 0.) then ! vapor - > ice
                tmp = tice - tz (k)
                ! 20160912: the following should produce more ice at higher altitude
                ! qi_crt = 4.92e-11 * exp (1.33 * log (1.e3 * exp (0.1 * tmp))) / den (k)
                qi_crt = qi_gen * min (qi_lim, 0.1 * tmp) / den (k)
                sink = min (sink, max (qi_crt - qi (k), pidep), tmp / tcpk (k))
            else ! ice -- > vapor
                pidep = pidep * min (1., dim (tz (k), t_sub) * 0.2)
                sink = max (pidep, sink, - qi (k))
            endif
            qv (k) = qv (k) - sink
            qi (k) = qi (k) + sink
            q_sol (k) = q_sol (k) + sink
            cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
            tz (k) = tz (k) + sink * (lhl (k) + lhi (k)) / cvm (k)
        endif

        ! -----------------------------------------------------------------------
        !> - Update capacity heat and latend heat coefficient.
        ! -----------------------------------------------------------------------

        lhl (k) = lv00 + d0_vap * tz (k)
        lhi (k) = li00 + dc_ice * tz (k)
        lcpk (k) = lhl (k) / cvm (k)
        icpk (k) = lhi (k) / cvm (k)
        tcpk (k) = lcpk (k) + icpk (k)

        ! -----------------------------------------------------------------------
        !> - Sublimation/deposition of snow.
        ! this process happens for all temp rage
        ! -----------------------------------------------------------------------

        if (qs (k) > qrmin) then
            qsi = iqs2 (tz (k), den (k), dqsdt)
            qden = qs (k) * den (k)
            tmp = exp (0.65625 * log (qden))
            tsq = tz (k) * tz (k)
            dq = (qsi - qv (k)) / (1. + tcpk (k) * dqsdt)
            pssub = cssub (1) * tsq * (cssub (2) * sqrt (qden) + cssub (3) * tmp * &
                sqrt (denfac (k))) / (cssub (4) * tsq + cssub (5) * qsi * den (k))
            pssub = (qsi - qv (k)) * dts * pssub
            if (pssub > 0.) then ! qs -- > qv, sublimation
                pssub = min (pssub * min (1., dim (tz (k), t_sub) * 0.2), qs (k))
            else
                if (tz (k) > tice) then
                    pssub = 0. ! no deposition
                else
                    pssub = max (pssub, dq, (tz (k) - tice) / tcpk (k))
                endif
            endif
            qs (k) = qs (k) - pssub
            qv (k) = qv (k) + pssub
            q_sol (k) = q_sol (k) - pssub
            cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
            tz (k) = tz (k) - pssub * (lhl (k) + lhi (k)) / cvm (k)
        endif

        ! -----------------------------------------------------------------------
        !> - Update capacity heat and latend heat coefficient.
        ! -----------------------------------------------------------------------

        lhl (k) = lv00 + d0_vap * tz (k)
        lhi (k) = li00 + dc_ice * tz (k)
        lcpk (k) = lhl (k) / cvm (k)
        icpk (k) = lhi (k) / cvm (k)
        tcpk (k) = lcpk (k) + icpk (k)

        ! -----------------------------------------------------------------------
        !> - Simplified 2-way grapuel sublimation-deposition mechanism.
        ! -----------------------------------------------------------------------

        if (qg (k) > qrmin) then
            qsi = iqs2 (tz (k), den (k), dqsdt)
            dq = (qv (k) - qsi) / (1. + tcpk (k) * dqsdt)
            pgsub = (qv (k) / qsi - 1.) * qg (k)
            if (pgsub > 0.) then ! deposition
                if (tz (k) > tice) then
                    pgsub = 0. ! no deposition
                else
                    pgsub = min (fac_v2g * pgsub, 0.2 * dq, ql (k) + qr (k), &
                        (tice - tz (k)) / tcpk (k))
                endif
            else ! submilation
                pgsub = max (fac_g2v * pgsub, dq) * min (1., dim (tz (k), t_sub) * 0.1)
            endif
            qg (k) = qg (k) + pgsub
            qv (k) = qv (k) - pgsub
            q_sol (k) = q_sol (k) + pgsub
            cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
            tz (k) = tz (k) + pgsub * (lhl (k) + lhi (k)) / cvm (k)
        endif

#ifdef USE_MIN_EVAP
        ! -----------------------------------------------------------------------
        !> - Update capacity heat and latend heat coefficient.
        ! -----------------------------------------------------------------------

        lhl (k) = lv00 + d0_vap * tz (k)
        lcpk (k) = lhl (k) / cvm (k)

        ! -----------------------------------------------------------------------
        !> - Minimum evap of rain in dry environmental air.
        ! -----------------------------------------------------------------------

        if (qr (k) > qcmin) then
            qsw = wqs2 (tz (k), den (k), dqsdt)
            sink = min (qr (k), dim (rh_rain * qsw, qv (k)) / (1. + lcpk (k) * dqsdt))
            qv (k) = qv (k) + sink
            qr (k) = qr (k) - sink
            q_liq (k) = q_liq (k) - sink
            cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
            tz (k) = tz (k) - sink * lhl (k) / cvm (k)
        endif
#endif

        ! -----------------------------------------------------------------------
        !> - Update capacity heat and latend heat coefficient.
        ! -----------------------------------------------------------------------

        lhl (k) = lv00 + d0_vap * tz (k)
        cvm (k) = c_air + (qv (k) + q_liq (k) + q_sol (k)) * c_vap
        lcpk (k) = lhl (k) / cvm (k)

        ! -----------------------------------------------------------------------
        ! compute cloud fraction
        ! -----------------------------------------------------------------------

        ! -----------------------------------------------------------------------
        !> - Combine water species.
        ! -----------------------------------------------------------------------

        if (do_qa) cycle

        if (rad_snow) then
            q_sol (k) = qi (k) + qs (k)
        else
            q_sol (k) = qi (k)
        endif
        if (rad_rain) then
            q_liq (k) = ql (k) + qr (k)
        else
            q_liq (k) = ql (k)
        endif
        q_cond (k) = q_liq (k) + q_sol (k)

        qpz = qv (k) + q_cond (k) ! qpz is conserved

        ! -----------------------------------------------------------------------
        !> - Use the "liquid-frozen water temperature" (tin) to compute saturated specific humidity.
        ! -----------------------------------------------------------------------

        tin = tz (k) - (lcpk (k) * q_cond (k) + icpk (k) * q_sol (k)) ! minimum temperature
        ! tin = tz (k) - ((lv00 + d0_vap * tz (k)) * q_cond (k) + &
        ! (li00 + dc_ice * tz (k)) * q_sol (k)) / (c_air + qpz * c_vap)

        ! -----------------------------------------------------------------------
        ! determine saturated specific humidity
        ! -----------------------------------------------------------------------

        if (tin <= t_wfr) then
            ! ice phase:
            qstar = iqs1 (tin, den (k))
        elseif (tin >= tice) then
            ! liquid phase:
            qstar = wqs1 (tin, den (k))
        else
            ! mixed phase:
            qsi = iqs1 (tin, den (k))
            qsw = wqs1 (tin, den (k))
            if (q_cond (k) > 3.e-6) then
                rqi = q_sol (k) / q_cond (k)
            else
                ! -----------------------------------------------------------------------
                ! mostly liquid water q_cond (k) at initial cloud development stage
                ! -----------------------------------------------------------------------
                rqi = (tice - tin) / (tice - t_wfr)
            endif
            qstar = rqi * qsi + (1. - rqi) * qsw
        endif

        ! -----------------------------------------------------------------------
        !> - Compute cloud fraction, assuming subgrid linear distribution in horizontal;
        !! this is effectively a smoother for the binary cloud scheme.
        ! -----------------------------------------------------------------------

        if (qpz > qrmin) then
            ! partial cloudiness by pdf:
            dq = max (qcmin, h_var * qpz)
            q_plus = qpz + dq ! cloud free if qstar > q_plus
            q_minus = qpz - dq
            if (qstar < q_minus) then
                qa (k) = qa (k) + 1. ! air fully saturated; 100 % cloud cover
            elseif (qstar < q_plus .and. q_cond (k) > qc_crt) then
                qa (k) = qa (k) + (q_plus - qstar) / (dq + dq) ! partial cloud cover
                ! qa (k) = sqrt (qa (k) + (q_plus - qstar) / (dq + dq))
            endif
        endif

    enddo

end subroutine subgrid_z_proc

! =======================================================================
!>\ingroup mod_gfdl_cloud_mp
!> This subroutine calculates rain evaporation.
subroutine revap_rac1 (hydrostatic, is, ie, dt, tz, qv, ql, qr, qi, qs, qg, den, hvar)

    implicit none

    logical, intent (in) :: hydrostatic

    integer, intent (in) :: is, ie

    real, intent (in) :: dt ! time step (s)

    real, intent (in), dimension (is:ie) :: den, hvar, qi, qs, qg

    real, intent (inout), dimension (is:ie) :: tz, qv, qr, ql

    real, dimension (is:ie) :: lcp2, denfac, q_liq, q_sol, cvm, lhl

    real :: dqv, qsat, dqsdt, evap, qden, q_plus, q_minus, sink
    real :: tin, t2, qpz, dq, dqh

    integer :: i

    ! -----------------------------------------------------------------------
    ! define latend heat coefficient
    ! -----------------------------------------------------------------------

    do i = is, ie
        lhl (i) = lv00 + d0_vap * tz (i)
        q_liq (i) = ql (i) + qr (i)
        q_sol (i) = qi (i) + qs (i) + qg (i)
        cvm (i) = c_air + qv (i) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice
        lcp2 (i) = lhl (i) / cvm (i)
        ! denfac (i) = sqrt (sfcrho / den (i))
    enddo

    do i = is, ie
        if (qr (i) > qrmin .and. tz (i) > t_wfr) then
            qpz = qv (i) + ql (i)
            tin = tz (i) - lcp2 (i) * ql (i) ! presence of clouds suppresses the rain evap
            qsat = wqs2 (tin, den (i), dqsdt)
            dqh = max (ql (i), hvar (i) * max (qpz, qcmin))
            dqv = qsat - qv (i)
            q_minus = qpz - dqh
            q_plus = qpz + dqh

            ! -----------------------------------------------------------------------
            ! qsat must be > q_minus to activate evaporation
            ! qsat must be < q_plus to activate accretion
            ! -----------------------------------------------------------------------

            ! -----------------------------------------------------------------------
            ! rain evaporation
            ! -----------------------------------------------------------------------

            if (dqv > qvmin .and. qsat > q_minus) then
                if (qsat > q_plus) then
                    dq = qsat - qpz
                else
                    ! q_minus < qsat < q_plus
                    ! dq == dqh if qsat == q_minus
                    dq = 0.25 * (q_minus - qsat) ** 2 / dqh
                endif
                qden = qr (i) * den (i)
                t2 = tin * tin
                evap = crevp (1) * t2 * dq * (crevp (2) * sqrt (qden) + crevp (3) * exp (0.725 * log (qden))) &
                     / (crevp (4) * t2 + crevp (5) * qsat * den (i))
                evap = min (qr (i), dt * evap, dqv / (1. + lcp2 (i) * dqsdt))
                qr (i) = qr (i) - evap
                qv (i) = qv (i) + evap
                q_liq (i) = q_liq (i) - evap
                cvm (i) = c_air + qv (i) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice
                tz (i) = tz (i) - evap * lhl (i) / cvm (i)
            endif

            ! -----------------------------------------------------------------------
            ! accretion: pracc
            ! -----------------------------------------------------------------------

            if (qr (i) > qrmin .and. ql (i) > 1.e-8 .and. qsat < q_plus) then
                denfac (i) = sqrt (sfcrho / den (i))
                sink = dt * denfac (i) * cracw * exp (0.95 * log (qr (i) * den (i)))
                sink = sink / (1. + sink) * ql (i)
                ql (i) = ql (i) - sink
                qr (i) = qr (i) + sink
            endif
        endif
    enddo

end subroutine revap_rac1

! =======================================================================
!>\ingroup mod_gfdl_cloud_mp
!>@brief The subroutine 'terminal_fall' computes terminal fall speed.
!>@details It considers cloud ice, snow, and graupel's melting during fall.
subroutine terminal_fall (dtm, ktop, kbot, tz, qv, ql, qr, qg, qs, qi, dz, dp, &
        den, vtg, vts, vti, r1, g1, s1, i1, m1_sol, w1)

    implicit none

    integer, intent (in) :: ktop, kbot

    real, intent (in) :: dtm ! time step (s)

    real, intent (in), dimension (ktop:kbot) :: vtg, vts, vti, den, dp, dz

    real, intent (inout), dimension (ktop:kbot) :: qv, ql, qr, qg, qs, qi, tz, m1_sol, w1

    real, intent (out) :: r1, g1, s1, i1

    real, dimension (ktop:kbot + 1) :: ze, zt

    real :: qsat, dqsdt, dt5, evap, dtime
    real :: factor, frac
    real :: tmp, precip, tc, sink

    real, dimension (ktop:kbot) :: lcpk, icpk, cvm, q_liq, q_sol, lhl, lhi
    real, dimension (ktop:kbot) :: m1, dm

    real :: zs = 0.
    real :: fac_imlt

    integer :: k, k0, m

    logical :: no_fall

    dt5 = 0.5 * dtm
    fac_imlt = 1. - exp (- dt5 / tau_imlt)

    ! -----------------------------------------------------------------------
    ! define heat capacity and latend heat coefficient
    ! -----------------------------------------------------------------------

    do k = ktop, kbot
        m1_sol (k) = 0.
        lhl (k) = lv00 + d0_vap * tz (k)
        lhi (k) = li00 + dc_ice * tz (k)
        q_liq (k) = ql (k) + qr (k)
        q_sol (k) = qi (k) + qs (k) + qg (k)
        cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
        lcpk (k) = lhl (k) / cvm (k)
        icpk (k) = lhi (k) / cvm (k)
    enddo

    ! -----------------------------------------------------------------------
    ! find significant melting level
    ! -----------------------------------------------------------------------

    k0 = kbot
    do k = ktop, kbot - 1
        if (tz (k) > tice) then
            k0 = k
            exit
        endif
    enddo

    ! -----------------------------------------------------------------------
    ! melting of cloud_ice (before fall) :
    ! -----------------------------------------------------------------------

    do k = k0, kbot
        tc = tz (k) - tice
        if (qi (k) > qcmin .and. tc > 0.) then
            sink = min (qi (k), fac_imlt * tc / icpk (k))
            tmp = min (sink, dim (ql_mlt, ql (k)))
            ql (k) = ql (k) + tmp
            qr (k) = qr (k) + sink - tmp
            qi (k) = qi (k) - sink
            q_liq (k) = q_liq (k) + sink
            q_sol (k) = q_sol (k) - sink
            cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
            tz (k) = tz (k) - sink * lhi (k) / cvm (k)
            tc = tz (k) - tice
        endif
    enddo

    ! -----------------------------------------------------------------------
    ! turn off melting when cloud microphysics time step is small
    ! -----------------------------------------------------------------------

    if (dtm < 60.) k0 = kbot

    ! sjl, turn off melting of falling cloud ice, snow and graupel
    k0 = kbot
    ! sjl, turn off melting of falling cloud ice, snow and graupel

    ze (kbot + 1) = zs
    do k = kbot, ktop, - 1
        ze (k) = ze (k + 1) - dz (k) ! dz < 0
    enddo

    zt (ktop) = ze (ktop)

    ! -----------------------------------------------------------------------
    ! update capacity heat and latend heat coefficient
    ! -----------------------------------------------------------------------

    do k = k0, kbot
        lhi (k) = li00 + dc_ice * tz (k)
        icpk (k) = lhi (k) / cvm (k)
    enddo

    ! -----------------------------------------------------------------------
    ! melting of falling cloud ice into rain
    ! -----------------------------------------------------------------------

    call check_column (ktop, kbot, qi, no_fall)

    if (vi_fac < 1.e-5 .or. no_fall) then
        i1 = 0.
    else

        do k = ktop + 1, kbot
            zt (k) = ze (k) - dt5 * (vti (k - 1) + vti (k))
        enddo
        zt (kbot + 1) = zs - dtm * vti (kbot)

        do k = ktop, kbot
            if (zt (k + 1) >= zt (k)) zt (k + 1) = zt (k) - dz_min
        enddo

        if (k0 < kbot) then
            do k = kbot - 1, k0, - 1
                if (qi (k) > qrmin) then
                    do m = k + 1, kbot
                        if (zt (k + 1) >= ze (m)) exit
                        if (zt (k) < ze (m + 1) .and. tz (m) > tice) then
                            dtime = min (1.0, (ze (m) - ze (m + 1)) / (max (vr_min, vti (k)) * tau_imlt))
                            sink = min (qi (k) * dp (k) / dp (m), dtime * (tz (m) - tice) / icpk (m))
                            tmp = min (sink, dim (ql_mlt, ql (m)))
                            ql (m) = ql (m) + tmp
                            qr (m) = qr (m) - tmp + sink
                            tz (m) = tz (m) - sink * icpk (m)
                            qi (k) = qi (k) - sink * dp (m) / dp (k)
                        endif
                    enddo
                endif
            enddo
        endif

        if (do_sedi_w) then
            do k = ktop, kbot
                dm (k) = dp (k) * (1. + qv (k) + ql (k) + qr (k) + qi (k) + qs (k) + qg (k))
            enddo
        endif

        if (use_ppm) then
            call lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, qi, i1, m1_sol, mono_prof)
        else
            call implicit_fall (dtm, ktop, kbot, ze, vti, dp, qi, i1, m1_sol)
        endif

        if (do_sedi_w) then
            w1 (ktop) = (dm (ktop) * w1 (ktop) + m1_sol (ktop) * vti (ktop)) / (dm (ktop) - m1_sol (ktop))
            do k = ktop + 1, kbot
                w1 (k) = (dm (k) * w1 (k) - m1_sol (k - 1) * vti (k - 1) + m1_sol (k) * vti (k)) &
                     / (dm (k) + m1_sol (k - 1) - m1_sol (k))
            enddo
        endif

    endif

    ! -----------------------------------------------------------------------
    ! melting of falling snow into rain
    ! -----------------------------------------------------------------------

    r1 = 0.

    call check_column (ktop, kbot, qs, no_fall)

    if (no_fall) then
        s1 = 0.
    else

        do k = ktop + 1, kbot
            zt (k) = ze (k) - dt5 * (vts (k - 1) + vts (k))
        enddo
        zt (kbot + 1) = zs - dtm * vts (kbot)

        do k = ktop, kbot
            if (zt (k + 1) >= zt (k)) zt (k + 1) = zt (k) - dz_min
        enddo

        if (k0 < kbot) then
            do k = kbot - 1, k0, - 1
                if (qs (k) > qrmin) then
                    do m = k + 1, kbot
                        if (zt (k + 1) >= ze (m)) exit
                        dtime = min (dtm, (ze (m) - ze (m + 1)) / (vr_min + vts (k)))
                        if (zt (k) < ze (m + 1) .and. tz (m) > tice) then
                            dtime = min (1.0, dtime / tau_smlt)
                            sink = min (qs (k) * dp (k) / dp (m), dtime * (tz (m) - tice) / icpk (m))
                            tz (m) = tz (m) - sink * icpk (m)
                            qs (k) = qs (k) - sink * dp (m) / dp (k)
                            if (zt (k) < zs) then
                                r1 = r1 + sink * dp (m) ! precip as rain
                            else
                                ! qr source here will fall next time step (therefore, can evap)
                                qr (m) = qr (m) + sink
                            endif
                        endif
                        if (qs (k) < qrmin) exit
                    enddo
                endif
            enddo
        endif

        if (do_sedi_w) then
            do k = ktop, kbot
                dm (k) = dp (k) * (1. + qv (k) + ql (k) + qr (k) + qi (k) + qs (k) + qg (k))
            enddo
        endif

        if (use_ppm) then
            call lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, qs, s1, m1, mono_prof)
        else
            call implicit_fall (dtm, ktop, kbot, ze, vts, dp, qs, s1, m1)
        endif

        do k = ktop, kbot
            m1_sol (k) = m1_sol (k) + m1 (k)
        enddo

        if (do_sedi_w) then
            w1 (ktop) = (dm (ktop) * w1 (ktop) + m1 (ktop) * vts (ktop)) / (dm (ktop) - m1 (ktop))
            do k = ktop + 1, kbot
                w1 (k) = (dm (k) * w1 (k) - m1 (k - 1) * vts (k - 1) + m1 (k) * vts (k)) &
                     / (dm (k) + m1 (k - 1) - m1 (k))
            enddo
        endif

    endif

    ! ----------------------------------------------
    ! melting of falling graupel into rain
    ! ----------------------------------------------

    call check_column (ktop, kbot, qg, no_fall)

    if (no_fall) then
        g1 = 0.
    else

        do k = ktop + 1, kbot
            zt (k) = ze (k) - dt5 * (vtg (k - 1) + vtg (k))
        enddo
        zt (kbot + 1) = zs - dtm * vtg (kbot)

        do k = ktop, kbot
            if (zt (k + 1) >= zt (k)) zt (k + 1) = zt (k) - dz_min
        enddo

        if (k0 < kbot) then
            do k = kbot - 1, k0, - 1
                if (qg (k) > qrmin) then
                    do m = k + 1, kbot
                        if (zt (k + 1) >= ze (m)) exit
                        dtime = min (dtm, (ze (m) - ze (m + 1)) / vtg (k))
                        if (zt (k) < ze (m + 1) .and. tz (m) > tice) then
                            dtime = min (1., dtime / tau_g2r)
                            sink = min (qg (k) * dp (k) / dp (m), dtime * (tz (m) - tice) / icpk (m))
                            tz (m) = tz (m) - sink * icpk (m)
                            qg (k) = qg (k) - sink * dp (m) / dp (k)
                            if (zt (k) < zs) then
                                r1 = r1 + sink * dp (m)
                            else
                                qr (m) = qr (m) + sink
                            endif
                        endif
                        if (qg (k) < qrmin) exit
                    enddo
                endif
            enddo
        endif

        if (do_sedi_w) then
            do k = ktop, kbot
                dm (k) = dp (k) * (1. + qv (k) + ql (k) + qr (k) + qi (k) + qs (k) + qg (k))
            enddo
        endif

        if (use_ppm) then
            call lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, qg, g1, m1, mono_prof)
        else
            call implicit_fall (dtm, ktop, kbot, ze, vtg, dp, qg, g1, m1)
        endif

        do k = ktop, kbot
            m1_sol (k) = m1_sol (k) + m1 (k)
        enddo

        if (do_sedi_w) then
            w1 (ktop) = (dm (ktop) * w1 (ktop) + m1 (ktop) * vtg (ktop)) / (dm (ktop) - m1 (ktop))
            do k = ktop + 1, kbot
                w1 (k) = (dm (k) * w1 (k) - m1 (k - 1) * vtg (k - 1) + m1 (k) * vtg (k)) &
                     / (dm (k) + m1 (k - 1) - m1 (k))
            enddo
        endif

    endif

end subroutine terminal_fall

! =======================================================================
!>\ingroup mod_gfdl_cloud_mp
!>@brief The subroutine 'check_column' checks
!! if the water species is large enough to fall.
subroutine check_column (ktop, kbot, q, no_fall)

    implicit none

    integer, intent (in) :: ktop, kbot

    real, intent (in) :: q (ktop:kbot)

    logical, intent (out) :: no_fall

    integer :: k

    no_fall = .true.

    do k = ktop, kbot
        if (q (k) > qrmin) then
            no_fall = .false.
            exit
        endif
    enddo

end subroutine check_column

! =======================================================================
!>\ingroup mod_gfdl_cloud_mp
!>@brief The subroutine computes the time-implicit monotonic
!! fall scheme.
!>@author Shian-Jiann Lin, 2016
subroutine implicit_fall (dt, ktop, kbot, ze, vt, dp, q, precip, m1)

    implicit none

    integer, intent (in) :: ktop, kbot

    real, intent (in) :: dt

    real, intent (in), dimension (ktop:kbot + 1) :: ze

    real, intent (in), dimension (ktop:kbot) :: vt, dp

    real, intent (inout), dimension (ktop:kbot) :: q

    real, intent (out), dimension (ktop:kbot) :: m1

    real, intent (out) :: precip

    real, dimension (ktop:kbot) :: dz, qm, dd

    integer :: k

    do k = ktop, kbot
        dz (k) = ze (k) - ze (k + 1)
        dd (k) = dt * vt (k)
        q (k) = q (k) * dp (k)
    enddo

    ! -----------------------------------------------------------------------
    ! sedimentation: non - vectorizable loop
    ! -----------------------------------------------------------------------

    qm (ktop) = q (ktop) / (dz (ktop) + dd (ktop))
    do k = ktop + 1, kbot
        qm (k) = (q (k) + dd (k - 1) * qm (k - 1)) / (dz (k) + dd (k))
    enddo

    ! -----------------------------------------------------------------------
    ! qm is density at this stage
    ! -----------------------------------------------------------------------

    do k = ktop, kbot
        qm (k) = qm (k) * dz (k)
    enddo

    ! -----------------------------------------------------------------------
    ! output mass fluxes: non - vectorizable loop
    ! -----------------------------------------------------------------------

    m1 (ktop) = q (ktop) - qm (ktop)
    do k = ktop + 1, kbot
        m1 (k) = m1 (k - 1) + q (k) - qm (k)
    enddo
    precip = m1 (kbot)

    ! -----------------------------------------------------------------------
    ! update:
    ! -----------------------------------------------------------------------

    do k = ktop, kbot
        q (k) = qm (k) / dp (k)
    enddo

end subroutine implicit_fall

! =======================================================================
!>\ingroup mod_gfdl_cloud_mp
!! Lagrangian scheme
!>  \author  S.J. Lin
subroutine lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, q, precip, m1, mono)

    implicit none

    integer, intent (in) :: ktop, kbot

    real, intent (in) :: zs

    logical, intent (in) :: mono

    real, intent (in), dimension (ktop:kbot + 1) :: ze, zt

    real, intent (in), dimension (ktop:kbot) :: dp

    ! m1: flux
    real, intent (inout), dimension (ktop:kbot) :: q, m1

    real, intent (out) :: precip

    real, dimension (ktop:kbot) :: qm, dz

    real :: a4 (4, ktop:kbot)

    real :: pl, pr, delz, esl

    integer :: k, k0, n, m

    real, parameter :: r3 = 1. / 3., r23 = 2. / 3.

    ! -----------------------------------------------------------------------
    ! density:
    ! -----------------------------------------------------------------------

    do k = ktop, kbot
        dz (k) = zt (k) - zt (k + 1) ! note: dz is positive
        q (k) = q (k) * dp (k)
        a4 (1, k) = q (k) / dz (k)
        qm (k) = 0.
    enddo

    ! -----------------------------------------------------------------------
    ! construct vertical profile with zt as coordinate
    ! -----------------------------------------------------------------------

    call cs_profile (a4 (1, ktop), dz (ktop), kbot - ktop + 1, mono)

    k0 = ktop
    do k = ktop, kbot
        do n = k0, kbot
            if (ze (k) <= zt (n) .and. ze (k) >= zt (n + 1)) then
                pl = (zt (n) - ze (k)) / dz (n)
                if (zt (n + 1) <= ze (k + 1)) then
                    ! entire new grid is within the original grid
                    pr = (zt (n) - ze (k + 1)) / dz (n)
                    qm (k) = a4 (2, n) + 0.5 * (a4 (4, n) + a4 (3, n) - a4 (2, n)) * (pr + pl) - &
                        a4 (4, n) * r3 * (pr * (pr + pl) + pl ** 2)
                    qm (k) = qm (k) * (ze (k) - ze (k + 1))
                    k0 = n
                    goto 555
                else
                    qm (k) = (ze (k) - zt (n + 1)) * (a4 (2, n) + 0.5 * (a4 (4, n) + &
                        a4 (3, n) - a4 (2, n)) * (1. + pl) - a4 (4, n) * (r3 * (1. + pl * (1. + pl))))
                    if (n < kbot) then
                        do m = n + 1, kbot
                            ! locate the bottom edge: ze (k + 1)
                            if (ze (k + 1) < zt (m + 1)) then
                                qm (k) = qm (k) + q (m)
                            else
                                delz = zt (m) - ze (k + 1)
                                esl = delz / dz (m)
                                qm (k) = qm (k) + delz * (a4 (2, m) + 0.5 * esl * &
                                     (a4 (3, m) - a4 (2, m) + a4 (4, m) * (1. - r23 * esl)))
                                k0 = m
                                goto 555
                            endif
                        enddo
                    endif
                    goto 555
                endif
            endif
        enddo
        555 continue
    enddo

    m1 (ktop) = q (ktop) - qm (ktop)
    do k = ktop + 1, kbot
        m1 (k) = m1 (k - 1) + q (k) - qm (k)
    enddo
    precip = m1 (kbot)

    ! convert back to * dry * mixing ratio:
    ! dp must be dry air_mass (because moist air mass will be changed due to terminal fall) .

    do k = ktop, kbot
        q (k) = qm (k) / dp (k)
    enddo

end subroutine lagrangian_fall_ppm

!>\ingroup mod_gfdl_cloud_mp
subroutine cs_profile (a4, del, km, do_mono)

    implicit none

    integer, intent (in) :: km ! vertical dimension

    real, intent (in) :: del (km)

    logical, intent (in) :: do_mono

    real, intent (inout) :: a4 (4, km)

    real, parameter :: qp_min = 1.e-6

    real :: gam (km)
    real :: q (km + 1)
    real :: d4, bet, a_bot, grat, pmp, lac
    real :: pmp_1, lac_1, pmp_2, lac_2
    real :: da1, da2, a6da

    integer :: k

    logical extm (km)

    grat = del (2) / del (1) ! grid ratio
    bet = grat * (grat + 0.5)
    q (1) = (2. * grat * (grat + 1.) * a4 (1, 1) + a4 (1, 2)) / bet
    gam (1) = (1. + grat * (grat + 1.5)) / bet

    do k = 2, km
        d4 = del (k - 1) / del (k)
        bet = 2. + 2. * d4 - gam (k - 1)
        q (k) = (3. * (a4 (1, k - 1) + d4 * a4 (1, k)) - q (k - 1)) / bet
        gam (k) = d4 / bet
    enddo

    a_bot = 1. + d4 * (d4 + 1.5)
    q (km + 1) = (2. * d4 * (d4 + 1.) * a4 (1, km) + a4 (1, km - 1) - a_bot * q (km)) &
         / (d4 * (d4 + 0.5) - a_bot * gam (km))

    do k = km, 1, - 1
        q (k) = q (k) - gam (k) * q (k + 1)
    enddo

    ! -----------------------------------------------------------------------
    ! apply constraints
    ! -----------------------------------------------------------------------

    do k = 2, km
        gam (k) = a4 (1, k) - a4 (1, k - 1)
    enddo

    ! -----------------------------------------------------------------------
    ! apply large - scale constraints to all fields if not local max / min
    ! -----------------------------------------------------------------------

    ! -----------------------------------------------------------------------
    ! top:
    ! -----------------------------------------------------------------------

    q (1) = max (q (1), 0.)
    q (2) = min (q (2), max (a4 (1, 1), a4 (1, 2)))
    q (2) = max (q (2), min (a4 (1, 1), a4 (1, 2)), 0.)

    ! -----------------------------------------------------------------------
    ! interior:
    ! -----------------------------------------------------------------------

    do k = 3, km - 1
        if (gam (k - 1) * gam (k + 1) > 0.) then
            q (k) = min (q (k), max (a4 (1, k - 1), a4 (1, k)))
            q (k) = max (q (k), min (a4 (1, k - 1), a4 (1, k)))
        else
            if (gam (k - 1) > 0.) then
                ! there exists a local max
                q (k) = max (q (k), min (a4 (1, k - 1), a4 (1, k)))
            else
                ! there exists a local min
                q (k) = min (q (k), max (a4 (1, k - 1), a4 (1, k)))
                q (k) = max (q (k), 0.0)
            endif
        endif
    enddo

    ! -----------------------------------------------------------------------
    ! bottom :
    ! -----------------------------------------------------------------------

    q (km) = min (q (km), max (a4 (1, km - 1), a4 (1, km)))
    q (km) = max (q (km), min (a4 (1, km - 1), a4 (1, km)), 0.)
    ! q (km + 1) = max (q (km + 1), 0.)

    ! -----------------------------------------------------------------------
    ! f (s) = al + s * [ (ar - al) + a6 * (1 - s) ] (0 <= s <= 1)
    ! -----------------------------------------------------------------------

    do k = 1, km - 1
        a4 (2, k) = q (k)
        a4 (3, k) = q (k + 1)
    enddo

    do k = 2, km - 1
        if (gam (k) * gam (k + 1) > 0.0) then
            extm (k) = .false.
        else
            extm (k) = .true.
        endif
    enddo

    if (do_mono) then
        do k = 3, km - 2
            if (extm (k)) then
                ! positive definite constraint only if true local extrema
                if (a4 (1, k) < qp_min .or. extm (k - 1) .or. extm (k + 1)) then
                    a4 (2, k) = a4 (1, k)
                    a4 (3, k) = a4 (1, k)
                endif
            else
                a4 (4, k) = 6. * a4 (1, k) - 3. * (a4 (2, k) + a4 (3, k))
                if (abs (a4 (4, k)) > abs (a4 (2, k) - a4 (3, k))) then
                    ! check within the smooth region if subgrid profile is non - monotonic
                    pmp_1 = a4 (1, k) - 2.0 * gam (k + 1)
                    lac_1 = pmp_1 + 1.5 * gam (k + 2)
                    a4 (2, k) = min (max (a4 (2, k), min (a4 (1, k), pmp_1, lac_1)), &
                        max (a4 (1, k), pmp_1, lac_1))
                    pmp_2 = a4 (1, k) + 2.0 * gam (k)
                    lac_2 = pmp_2 - 1.5 * gam (k - 1)
                    a4 (3, k) = min (max (a4 (3, k), min (a4 (1, k), pmp_2, lac_2)), &
                        max (a4 (1, k), pmp_2, lac_2))
                endif
            endif
        enddo
    else
        do k = 3, km - 2
            if (extm (k)) then
                if (a4 (1, k) < qp_min .or. extm (k - 1) .or. extm (k + 1)) then
                    a4 (2, k) = a4 (1, k)
                    a4 (3, k) = a4 (1, k)
                endif
            endif
        enddo
    endif

    do k = 1, km - 1
        a4 (4, k) = 6. * a4 (1, k) - 3. * (a4 (2, k) + a4 (3, k))
    enddo

    k = km - 1
    if (extm (k)) then
        a4 (2, k) = a4 (1, k)
        a4 (3, k) = a4 (1, k)
        a4 (4, k) = 0.
    else
        da1 = a4 (3, k) - a4 (2, k)
        da2 = da1 ** 2
        a6da = a4 (4, k) * da1
        if (a6da < - da2) then
            a4 (4, k) = 3. * (a4 (2, k) - a4 (1, k))
            a4 (3, k) = a4 (2, k) - a4 (4, k)
        elseif (a6da > da2) then
            a4 (4, k) = 3. * (a4 (3, k) - a4 (1, k))
            a4 (2, k) = a4 (3, k) - a4 (4, k)
        endif
    endif

    call cs_limiters (km - 1, a4)

    ! -----------------------------------------------------------------------
    ! bottom layer:
    ! -----------------------------------------------------------------------

    a4 (2, km) = a4 (1, km)
    a4 (3, km) = a4 (1, km)
    a4 (4, km) = 0.

end subroutine cs_profile

!>\ingroup mod_gfdl_cloud_mp
!! This subroutine perform positive definite constraint.
subroutine cs_limiters (km, a4)

    implicit none

    integer, intent (in) :: km

    real, intent (inout) :: a4 (4, km) ! ppm array

    real, parameter :: r12 = 1. / 12.

    integer :: k

    ! -----------------------------------------------------------------------
    ! positive definite constraint
    ! -----------------------------------------------------------------------

    do k = 1, km
        if (abs (a4 (3, k) - a4 (2, k)) < - a4 (4, k)) then
            if ((a4 (1, k) + 0.25 * (a4 (3, k) - a4 (2, k)) ** 2 / a4 (4, k) + a4 (4, k) * r12) < 0.) then
                if (a4 (1, k) < a4 (3, k) .and. a4 (1, k) < a4 (2, k)) then
                    a4 (3, k) = a4 (1, k)
                    a4 (2, k) = a4 (1, k)
                    a4 (4, k) = 0.
                elseif (a4 (3, k) > a4 (2, k)) then
                    a4 (4, k) = 3. * (a4 (2, k) - a4 (1, k))
                    a4 (3, k) = a4 (2, k) - a4 (4, k)
                else
                    a4 (4, k) = 3. * (a4 (3, k) - a4 (1, k))
                    a4 (2, k) = a4 (3, k) - a4 (4, k)
                endif
            endif
        endif
    enddo

end subroutine cs_limiters

! =======================================================================
!>\ingroup mod_gfdl_cloud_mp
!! The subroutine calculates vertical fall speed of snow/ice/graupel.
subroutine fall_speed (ktop, kbot, den, qs, qi, qg, ql, tk, vts, vti, vtg)

    implicit none

    integer, intent (in) :: ktop, kbot

    real, intent (in), dimension (ktop:kbot) :: den, qs, qi, qg, ql, tk
    real, intent (out), dimension (ktop:kbot) :: vts, vti, vtg

    ! fall velocity constants:

    real, parameter :: thi = 1.0e-8 !< cloud ice threshold for terminal fall
    real, parameter :: thg = 1.0e-8
    real, parameter :: ths = 1.0e-8

    ! coefficient for the parameterization of mass weighted fall velocity
    ! as a function of temperature and IWC.
    ! Table 1 in Deng and Mace (2008) \cite deng_and_mace_2008.
    real, parameter :: aa = - 4.14122e-5
    real, parameter :: bb = - 0.00538922
    real, parameter :: cc = - 0.0516344
    real, parameter :: dd = 0.00216078
    real, parameter :: ee = 1.9714

    ! marshall - palmer constants

    real, parameter :: vcons = 6.6280504
    real, parameter :: vcong = 87.2382675
    real, parameter :: vconh = vcong * sqrt (rhoh / rhog)
    real, parameter :: norms = 942477796.076938
    real, parameter :: normg = 5026548245.74367
    real, parameter :: normh = pi * rhoh * rnzh

    real, dimension (ktop:kbot) :: qden, tc, rhof

    real :: vi0

    integer :: k

    ! -----------------------------------------------------------------------
    ! marshall - palmer formula
    ! -----------------------------------------------------------------------

    ! -----------------------------------------------------------------------
    ! try the local air density -- for global model; the true value could be
    ! much smaller than sfcrho over high mountains
    ! -----------------------------------------------------------------------

    do k = ktop, kbot
        rhof (k) = sqrt (min (10., sfcrho / den (k)))
    enddo

    ! -----------------------------------------------------------------------
    !> - ice: use Deng and Mace (2008) \cite deng_and_mace_2008, which gives smaler
    !! fall speed than Heymsfield and Donner (1990) \cite heymsfield_and_donner_1990.
    ! -----------------------------------------------------------------------

    if (const_vi) then
        vti (:) = vi_fac
    else
        ! -----------------------------------------------------------------------
        ! use deng and mace (2008, grl), which gives smaller fall speed than hd90 formula
        ! -----------------------------------------------------------------------
        vi0 = 0.01 * vi_fac
        do k = ktop, kbot
            if (qi (k) < thi) then ! this is needed as the fall - speed maybe problematic for small qi
                vti (k) = vf_min
            else
                tc (k) = tk (k) - tice
                vti (k) = (3. + log10 (qi (k) * den (k))) * (tc (k) * (aa * tc (k) + bb) + cc) + dd * tc (k) + ee
                vti (k) = vi0 * exp (log_10 * vti (k)) * 0.9
                vti (k) = min (vi_max, max (vf_min, vti (k)))
            endif
        enddo
    endif

    ! -----------------------------------------------------------------------
    !> - snow:
    ! -----------------------------------------------------------------------

    if (const_vs) then
        vts (:) = vs_fac ! 1. ifs_2016
    else
        do k = ktop, kbot
            if (qs (k) < ths) then
                vts (k) = vf_min
            else
                vts (k) = vs_fac * vcons * rhof (k) * exp (0.0625 * log (qs (k) * den (k) / norms))
                vts (k) = min (vs_max, max (vf_min, vts (k)))
            endif
        enddo
    endif

    ! -----------------------------------------------------------------------
    !> - graupel:
    ! -----------------------------------------------------------------------
    if (const_vg) then
        vtg (:) = vg_fac ! 2.
    else
        if (do_hail) then
            do k = ktop, kbot
                if (qg (k) < thg) then
                    vtg (k) = vf_min
                else
                    vtg (k) = vg_fac * vconh * rhof (k) * sqrt (sqrt (sqrt (qg (k) * den (k) / normh)))
                    vtg (k) = min (vg_max, max (vf_min, vtg (k)))
                endif
            enddo
        else
        do k = ktop, kbot
            if (qg (k) < thg) then
                vtg (k) = vf_min
            else
                vtg (k) = vg_fac * vcong * rhof (k) * sqrt (sqrt (sqrt (qg (k) * den (k) / normg)))
                vtg (k) = min (vg_max, max (vf_min, vtg (k)))
            endif
        enddo
    endif
    endif

end subroutine fall_speed

! =======================================================================
!>\ingroup mod_gfdl_cloud_mp
!! The subroutine sets up gfdl cloud microphysics parameters.
subroutine setupm

    implicit none

    real :: gcon, cd, scm3, pisq, act (8)
    real :: vdifu, tcond
    real :: visk
    real :: ch2o, hltf
    real :: hlts, hltc, ri50

    real, parameter :: gam263 = 1.456943, gam275 = 1.608355, gam290 = 1.827363, &
        gam325 = 2.54925, gam350 = 3.323363, gam380 = 4.694155, &
        gam425 = 8.285063, gam450 = 11.631769, gam480 = 17.837789, &
        gam625 = 184.860962, gam680 = 496.604067

    ! intercept parameters

!    real, parameter :: rnzr = 8.0e6 ! lin83
!    real, parameter :: rnzs = 3.0e6 ! lin83
!    real, parameter :: rnzg = 4.0e6 ! rh84

    ! density parameters

!    real, parameter :: rhos = 0.1e3 !< lin83 (snow density; 1 / 10 of water)
!    real, parameter :: rhog = 0.4e3 !< rh84 (graupel density)
    real, parameter :: acc (3) = (/ 5.0, 2.0, 0.5 /)

    real den_rc

    integer :: i, k

    pie = 4. * atan (1.0)

    ! s. klein's formular (eq 16) from am2

    fac_rc = (4. / 3.) * pie * rhor * rthresh ** 3

    if (prog_ccn) then
        ! if (master) write (*, *) 'prog_ccn option is .t.'
    else
        den_rc = fac_rc * ccn_o * 1.e6
        ! if (master) write (*, *) 'mp: for ccn_o = ', ccn_o, 'ql_rc = ', den_rc
        den_rc = fac_rc * ccn_l * 1.e6
        ! if (master) write (*, *) 'mp: for ccn_l = ', ccn_l, 'ql_rc = ', den_rc
    endif

    vdifu = 2.11e-5
    tcond = 2.36e-2

    visk = 1.259e-5
    hlts = 2.8336e6
    hltc = 2.5e6
    hltf = 3.336e5

    ch2o = 4.1855e3
    ri50 = 1.e-4

    pisq = pie * pie
    scm3 = (visk / vdifu) ** (1. / 3.)

    cracs = pisq * rnzr * rnzs * rhos
    csacr = pisq * rnzr * rnzs * rhor
    if (do_hail) then
        cgacr = pisq * rnzr * rnzh * rhor
        cgacs = pisq * rnzh * rnzs * rhos
    else
        cgacr = pisq * rnzr * rnzg * rhor
        cgacs = pisq * rnzg * rnzs * rhos
    endif
    cgacs = cgacs * c_pgacs

    ! act: 1 - 2:racs (s - r) ; 3 - 4:sacr (r - s) ;
    ! 5 - 6:gacr (r - g) ; 7 - 8:gacs (s - g)

    act (1) = pie * rnzs * rhos
    act (2) = pie * rnzr * rhor
    if (do_hail) then
        act (6) = pie * rnzh * rhoh
    else
        act (6) = pie * rnzg * rhog
    endif
    act (3) = act (2)
    act (4) = act (1)
    act (5) = act (2)
    act (7) = act (1)
    act (8) = act (6)

    do i = 1, 3
        do k = 1, 4
            acco (i, k) = acc (i) / (act (2 * k - 1) ** ((7 - i) * 0.25) * act (2 * k) ** (i * 0.25))
        enddo
    enddo

    gcon = 40.74 * sqrt (sfcrho) ! 44.628

    csacw = pie * rnzs * clin * gam325 / (4. * act (1) ** 0.8125)
    ! decreasing csacw to reduce cloud water --- > snow

    craci = pie * rnzr * alin * gam380 / (4. * act (2) ** 0.95)
    csaci = csacw * c_psaci

    if (do_hail) then
        cgacw = pie * rnzh * gam350 * gcon / (4. * act (6) ** 0.875)
    else
        cgacw = pie * rnzg * gam350 * gcon / (4. * act (6) ** 0.875)
    endif
    ! cgaci = cgacw * 0.1

    ! sjl, may 28, 2012
    cgaci = cgacw * 0.05
    ! sjl, may 28, 2012

    cracw = craci ! cracw = 3.27206196043822
    cracw = c_cracw * cracw

    ! subl and revp: five constants for three separate processes

    cssub (1) = 2. * pie * vdifu * tcond * rvgas * rnzs
    if (do_hail) then
        cgsub (1) = 2. * pie * vdifu * tcond * rvgas * rnzh
    else
        cgsub (1) = 2. * pie * vdifu * tcond * rvgas * rnzg
    endif
    crevp (1) = 2. * pie * vdifu * tcond * rvgas * rnzr
    cssub (2) = 0.78 / sqrt (act (1))
    cgsub (2) = 0.78 / sqrt (act (6))
    crevp (2) = 0.78 / sqrt (act (2))
    cssub (3) = 0.31 * scm3 * gam263 * sqrt (clin / visk) / act (1) ** 0.65625
    cgsub (3) = 0.31 * scm3 * gam275 * sqrt (gcon / visk) / act (6) ** 0.6875
    crevp (3) = 0.31 * scm3 * gam290 * sqrt (alin / visk) / act (2) ** 0.725
    cssub (4) = tcond * rvgas
    cssub (5) = hlts ** 2 * vdifu
    cgsub (4) = cssub (4)
    crevp (4) = cssub (4)
    cgsub (5) = cssub (5)
    crevp (5) = hltc ** 2 * vdifu

    cgfr (1) = 20.e2 * pisq * rnzr * rhor / act (2) ** 1.75
    cgfr (2) = 0.66

    ! smlt: five constants (lin et al. 1983)

    csmlt (1) = 2. * pie * tcond * rnzs / hltf
    csmlt (2) = 2. * pie * vdifu * rnzs * hltc / hltf
    csmlt (3) = cssub (2)
    csmlt (4) = cssub (3)
    csmlt (5) = ch2o / hltf

    ! gmlt: five constants

    if (do_hail) then
        cgmlt (1) = 2. * pie * tcond * rnzh / hltf
        cgmlt (2) = 2. * pie * vdifu * rnzh * hltc / hltf
    else
        cgmlt (1) = 2. * pie * tcond * rnzg / hltf
        cgmlt (2) = 2. * pie * vdifu * rnzg * hltc / hltf
    endif
    cgmlt (3) = cgsub (2)
    cgmlt (4) = cgsub (3)
    cgmlt (5) = ch2o / hltf

    es0 = 6.107799961e2 ! ~6.1 mb
    ces0 = eps * es0

end subroutine setupm

! =======================================================================
! initialization of gfdl cloud microphysics
!>\ingroup mod_gfdl_cloud_mp
!! The subroutine 'gfdl_cloud_microphys_init' initializes the GFDL
!! cloud microphysics.
subroutine gfdl_cloud_microphys_mod_init (me, master, nlunit, input_nml_file, logunit, fn_nml)

    implicit none

    integer, intent (in) :: me
    integer, intent (in) :: master
    integer, intent (in) :: nlunit
    integer, intent (in) :: logunit

    character (len = 64), intent (in) :: fn_nml
    character (len = *),  intent (in) :: input_nml_file(:)

    integer :: ios
    logical :: exists

    ! integer, intent (in) :: id, jd, kd
    ! integer, intent (in) :: axes (4)
    ! type (time_type), intent (in) :: time

    ! integer :: unit, io, ierr, k, logunit
    ! logical :: flag
    ! real :: tmp, q1, q2

    ! master = (mpp_pe () .eq.mpp_root_pe ())

#ifdef INTERNAL_FILE_NML
    read (input_nml_file, nml = gfdl_cloud_microphysics_nml)
#else
    inquire (file = trim (fn_nml), exist = exists)
    if (.not. exists) then
        write (6, *) 'gfdl - mp :: namelist file: ', trim (fn_nml), ' does not exist'
        stop
    else
        open (unit = nlunit, file = fn_nml, action = 'read' , status = 'old', iostat = ios)
    endif
    rewind (nlunit)
    read (nlunit, nml = gfdl_cloud_microphysics_nml)
    close (nlunit)
#endif

    ! write version number and namelist to log file
    if (me == master) then
        write (logunit, *) " ================================================================== "
        write (logunit, *) "gfdl_cloud_microphys_mod"
        write (logunit, nml = gfdl_cloud_microphysics_nml)
    endif

    if (do_setup) then
        call setup_con
        call setupm
        do_setup = .false.
    endif

    log_10 = log (10.)

    tice0 = tice - 0.01
    t_wfr = tice - 40.0 ! supercooled water can exist down to - 48 c, which is the "absolute"

    ! if (master) write (logunit, nml = gfdl_cloud_microphys_nml)
    !
    ! id_vtr = register_diag_field (mod_name, 'vt_r', axes (1:3), time, &
    ! 'rain fall speed', 'm / s', missing_value = missing_value)
    ! id_vts = register_diag_field (mod_name, 'vt_s', axes (1:3), time, &
    ! 'snow fall speed', 'm / s', missing_value = missing_value)
    ! id_vtg = register_diag_field (mod_name, 'vt_g', axes (1:3), time, &
    ! 'graupel fall speed', 'm / s', missing_value = missing_value)
    ! id_vti = register_diag_field (mod_name, 'vt_i', axes (1:3), time, &
    ! 'ice fall speed', 'm / s', missing_value = missing_value)

    ! id_droplets = register_diag_field (mod_name, 'droplets', axes (1:3), time, &
    ! 'droplet number concentration', '# / m3', missing_value = missing_value)
    ! id_rh = register_diag_field (mod_name, 'rh_lin', axes (1:2), time, &
    ! 'relative humidity', 'n / a', missing_value = missing_value)

    ! id_rain = register_diag_field (mod_name, 'rain_lin', axes (1:2), time, &
    ! 'rain_lin', 'mm / day', missing_value = missing_value)
    ! id_snow = register_diag_field (mod_name, 'snow_lin', axes (1:2), time, &
    ! 'snow_lin', 'mm / day', missing_value = missing_value)
    ! id_graupel = register_diag_field (mod_name, 'graupel_lin', axes (1:2), time, &
    ! 'graupel_lin', 'mm / day', missing_value = missing_value)
    ! id_ice = register_diag_field (mod_name, 'ice_lin', axes (1:2), time, &
    ! 'ice_lin', 'mm / day', missing_value = missing_value)
    ! id_prec = register_diag_field (mod_name, 'prec_lin', axes (1:2), time, &
    ! 'prec_lin', 'mm / day', missing_value = missing_value)

    ! if (master) write (*, *) 'prec_lin diagnostics initialized.', id_prec

    ! id_cond = register_diag_field (mod_name, 'cond_lin', axes (1:2), time, &
    ! 'total condensate', 'kg / m ** 2', missing_value = missing_value)
    ! id_var = register_diag_field (mod_name, 'var_lin', axes (1:2), time, &
    ! 'subgrid variance', 'n / a', missing_value = missing_value)

    ! call qsmith_init

    ! testing the water vapor tables

    ! if (mp_debug .and. master) then
    ! write (*, *) 'testing water vapor tables in gfdl_cloud_microphys'
    ! tmp = tice - 90.
    ! do k = 1, 25
    ! q1 = wqsat_moist (tmp, 0., 1.e5)
    ! q2 = qs1d_m (tmp, 0., 1.e5)
    ! write (*, *) nint (tmp - tice), q1, q2, 'dq = ', q1 - q2
    ! tmp = tmp + 5.
    ! enddo
    ! endif

    ! if (master) write (*, *) 'gfdl_cloud_micrphys diagnostics initialized.'

    module_is_initialized = .true.

!+---+-----------------------------------------------------------------+
!..Set these variables needed for computing radar reflectivity.  These
!.. get used within radar_init to create other variables used in the
!.. radar module.

   xam_r = pi*rhor/6.
   xbm_r = 3.
   xmu_r = 0.
   xam_s = pi*rhos/6.
   xbm_s = 3.
   xmu_s = 0.
   xam_g = pi*rhog/6.
   xbm_g = 3.
   xmu_g = 0.

   call radar_init

end subroutine gfdl_cloud_microphys_mod_init

! =======================================================================
! end of gfdl cloud microphysics
!>\ingroup mod_gfdl_cloud_mp
!! The subroutine 'gfdl_cloud_microphys_init' terminates the GFDL
!! cloud microphysics.
subroutine gfdl_cloud_microphys_mod_end()

    implicit none

    deallocate (table)
    deallocate (table2)
    deallocate (table3)
    deallocate (tablew)
    deallocate (des)
    deallocate (des2)
    deallocate (des3)
    deallocate (desw)

    tables_are_initialized = .false.

end subroutine gfdl_cloud_microphys_mod_end

! =======================================================================
! qsmith table initialization
!>\ingroup mod_gfdl_cloud_mp
!! The subroutine 'setup_con' sets up constants and calls 'qsmith_init'.
subroutine setup_con

    implicit none

    ! master = (mpp_pe () .eq.mpp_root_pe ())

    rgrav = 1. / grav

    if (.not. qsmith_tables_initialized) call qsmith_init

    qsmith_tables_initialized = .true.

end subroutine setup_con

! =======================================================================
!>\ingroup gfdlmp
!>@brief The function is an accretion function (Lin et al.(1983)
!! \cite lin_et_al_1983 )
! =======================================================================

real function acr3d (v1, v2, q1, q2, c, cac, rho)

    implicit none

    real, intent (in) :: v1, v2, c, rho
    real, intent (in) :: q1, q2 ! mixing ratio!!!
    real, intent (in) :: cac (3)

    real :: t1, s1, s2

    ! integer :: k
    !
    ! real :: a
    !
    ! a = 0.0
    ! do k = 1, 3
    ! a = a + cac (k) * ((q1 * rho) ** ((7 - k) * 0.25) * (q2 * rho) ** (k * 0.25))
    ! enddo
    ! acr3d = c * abs (v1 - v2) * a / rho

    ! optimized

    t1 = sqrt (q1 * rho)
    s1 = sqrt (q2 * rho)
    s2 = sqrt (s1) ! s1 = s2 ** 2
    acr3d = c * abs (v1 - v2) * q1 * s2 * (cac (1) * t1 + cac (2) * sqrt (t1) * s2 + cac (3) * s1)

end function acr3d

! =======================================================================
!>\ingroup gfdlmp
!>@brief Melting of snow function (Lin et al.(1983) \cite lin_et_al_1983)
!!  note: psacw and psacr must be calc before smlt is called
! =======================================================================

real function smlt (tc, dqs, qsrho, psacw, psacr, c, rho, rhofac)

    implicit none

    real, intent (in) :: tc, dqs, qsrho, psacw, psacr, c (5), rho, rhofac

    smlt = (c (1) * tc / rho - c (2) * dqs) * (c (3) * sqrt (qsrho) + &
        c (4) * qsrho ** 0.65625 * sqrt (rhofac)) + c (5) * tc * (psacw + psacr)

end function smlt

! =======================================================================
!>@brief Melting of graupel function (Eq.(47) in Lin et al. 1983 \cite lin_et_al_1983)
!!\n  note: \f$P_{gacw}\f$ and \f$P_{gacr}\f$ must be calculated before gmlt is called.
! =======================================================================

real function gmlt (tc, dqs, qgrho, pgacw, pgacr, c, rho)

    implicit none

    real, intent (in) :: tc, dqs, qgrho, pgacw, pgacr, c (5), rho

    gmlt = (c (1) * tc / rho - c (2) * dqs) * (c (3) * sqrt (qgrho) + &
        c (4) * qgrho ** 0.6875 / rho ** 0.25) + c (5) * tc * (pgacw + pgacr)

end function gmlt

! =======================================================================
! initialization
! prepare saturation water vapor pressure tables
! =======================================================================
!>\ingroup mod_gfdl_cloud_mp
!>@brief The subroutine 'qsmith_init' initializes lookup tables for saturation
!! water vapor pressure for the following utility routines that are designed
!! to return qs consistent with the assumptions in FV3.
!>@details The calculations are highly accurate values based on the Clausius-Clapeyron
!! equation.
! =======================================================================
subroutine qsmith_init

    implicit none

    integer, parameter :: length = 2621

    integer :: i

    if (.not. tables_are_initialized) then

       ! master = (mpp_pe () .eq. mpp_root_pe ())
       ! if (master) print *, ' gfdl mp: initializing qs tables'

       ! debug code
       ! print *, mpp_pe (), allocated (table), allocated (table2), &
       ! allocated (table3), allocated (tablew), allocated (des), &
       ! allocated (des2), allocated (des3), allocated (desw)
       ! end debug code

       ! generate es table (dt = 0.1 deg. c)

       allocate (table (length))
       allocate (table2 (length))
       allocate (table3 (length))
       allocate (tablew (length))
       allocate (des (length))
       allocate (des2 (length))
       allocate (des3 (length))
       allocate (desw (length))

       call qs_table (length)
       call qs_table2 (length)
       call qs_table3 (length)
       call qs_tablew (length)

       do i = 1, length - 1
           des (i) = max (0., table (i + 1) - table (i))
           des2 (i) = max (0., table2 (i + 1) - table2 (i))
           des3 (i) = max (0., table3 (i + 1) - table3 (i))
           desw (i) = max (0., tablew (i + 1) - tablew (i))
       enddo
       des (length) = des (length - 1)
       des2 (length) = des2 (length - 1)
       des3 (length) = des3 (length - 1)
       desw (length) = desw (length - 1)

        tables_are_initialized = .true.

    endif

end subroutine qsmith_init

! =======================================================================
! compute the saturated specific humidity for table ii
!>\ingroup mod_gfdl_cloud_mp
!>@brief The function 'wqs1' returns the saturation vapor pressure over pure
!! liquid water for a given temperature and air density.
real function wqs1 (ta, den)

    implicit none

    !> pure water phase; universal dry / moist formular using air density
    !> input "den" can be either dry or moist air density

    real, intent (in) :: ta, den

    real :: es, ap1, tmin

    integer :: it

    tmin = table_ice - 160.
    ap1 = 10. * dim (ta, tmin) + 1.
    ap1 = min (2621., ap1)
    it = ap1
    es = tablew (it) + (ap1 - it) * desw (it)
    wqs1 = es / (rvgas * ta * den)

end function wqs1

! =======================================================================
! compute the gradient of saturated specific humidity for table ii
!>@brief The function 'wqs2' returns the saturation vapor pressure over pure
!! liquid water for a given temperature and air density, as well as the
!! analytic dqs/dT: rate of change of saturation vapor pressure WRT temperature.
! =======================================================================

real function wqs2 (ta, den, dqdt)

    implicit none

    ! pure water phase; universal dry / moist formular using air density
    ! input "den" can be either dry or moist air density

    real, intent (in) :: ta, den

    real, intent (out) :: dqdt

    real :: es, ap1, tmin

    integer :: it

    tmin = table_ice - 160.

    if (.not. tables_are_initialized) call qsmith_init

    ap1 = 10. * dim (ta, tmin) + 1.
    ap1 = min (2621., ap1)
    it = ap1
    es = tablew (it) + (ap1 - it) * desw (it)
    wqs2 = es / (rvgas * ta * den)
    it = ap1 - 0.5
    ! finite diff, del_t = 0.1:
    dqdt = 10. * (desw (it) + (ap1 - it) * (desw (it + 1) - desw (it))) / (rvgas * ta * den)

end function wqs2

! =======================================================================
! compute wet buld temperature
!>@brief The function 'wet_bulb' uses 'wqs2' to compute the wet-bulb temperature
!! from the mixing ratio and the temperature.
! =======================================================================

real function wet_bulb (q, t, den)

    implicit none

    real, intent (in) :: t, q, den

    real :: qs, tp, dqdt

    wet_bulb = t
    qs = wqs2 (wet_bulb, den, dqdt)
    tp = 0.5 * (qs - q) / (1. + lcp * dqdt) * lcp
    wet_bulb = wet_bulb - tp

    ! tp is negative if super - saturated
    if (tp > 0.01) then
        qs = wqs2 (wet_bulb, den, dqdt)
        tp = (qs - q) / (1. + lcp * dqdt) * lcp
        wet_bulb = wet_bulb - tp
    endif

end function wet_bulb

! =======================================================================
!>@brief The function 'iqs1' computes the saturated specific humidity
!! for table iii
! =======================================================================

real function iqs1 (ta, den)

    implicit none

    !> water - ice phase; universal dry / moist formular using air density
    !> input "den" can be either dry or moist air density

    real, intent (in) :: ta, den

    real :: es, ap1, tmin

    integer :: it

    tmin = table_ice - 160.
    ap1 = 10. * dim (ta, tmin) + 1.
    ap1 = min (2621., ap1)
    it = ap1
    es = table2 (it) + (ap1 - it) * des2 (it)
    iqs1 = es / (rvgas * ta * den)

end function iqs1

! =======================================================================
!>@brief The function 'iqs2' computes the gradient of saturated specific
!! humidity for table iii
! =======================================================================

real function iqs2 (ta, den, dqdt)

    implicit none

    ! water - ice phase; universal dry / moist formular using air density
    ! input "den" can be either dry or moist air density

    real, intent (in) :: ta, den

    real, intent (out) :: dqdt

    real :: es, ap1, tmin

    integer :: it

    tmin = table_ice - 160.
    ap1 = 10. * dim (ta, tmin) + 1.
    ap1 = min (2621., ap1)
    it = ap1
    es = table2 (it) + (ap1 - it) * des2 (it)
    iqs2 = es / (rvgas * ta * den)
    it = ap1 - 0.5
    dqdt = 10. * (des2 (it) + (ap1 - it) * (des2 (it + 1) - des2 (it))) / (rvgas * ta * den)

end function iqs2

! =======================================================================
!>@brief The function 'qs1d_moist' computes the gradient of saturated
!! specific humidity for table iii.
! =======================================================================

real function qs1d_moist (ta, qv, pa, dqdt)

    implicit none

    real, intent (in) :: ta, pa, qv

    real, intent (out) :: dqdt

    real :: es, ap1, tmin, eps10

    integer :: it

    tmin = table_ice - 160.
    eps10 = 10. * eps
    ap1 = 10. * dim (ta, tmin) + 1.
    ap1 = min (2621., ap1)
    it = ap1
    es = table2 (it) + (ap1 - it) * des2 (it)
    qs1d_moist = eps * es * (1. + zvir * qv) / pa
    it = ap1 - 0.5
    dqdt = eps10 * (des2 (it) + (ap1 - it) * (des2 (it + 1) - des2 (it))) * (1. + zvir * qv) / pa

end function qs1d_moist

! =======================================================================
! compute the gradient of saturated specific humidity for table ii
!>@brief The function 'wqsat2_moist' computes the saturated specific humidity
!! for pure liquid water , as well as des/dT.
! =======================================================================

real function wqsat2_moist (ta, qv, pa, dqdt)

    implicit none

    real, intent (in) :: ta, pa, qv

    real, intent (out) :: dqdt

    real :: es, ap1, tmin, eps10

    integer :: it

    tmin = table_ice - 160.
    eps10 = 10. * eps
    ap1 = 10. * dim (ta, tmin) + 1.
    ap1 = min (2621., ap1)
    it = ap1
    es = tablew (it) + (ap1 - it) * desw (it)
    wqsat2_moist = eps * es * (1. + zvir * qv) / pa
    it = ap1 - 0.5
    dqdt = eps10 * (desw (it) + (ap1 - it) * (desw (it + 1) - desw (it))) * (1. + zvir * qv) / pa

end function wqsat2_moist

! =======================================================================
! compute the saturated specific humidity for table ii
!>@brief The function 'wqsat_moist' computes the saturated specific humidity
!! for pure liquid water.
! =======================================================================

real function wqsat_moist (ta, qv, pa)

    implicit none

    real, intent (in) :: ta, pa, qv

    real :: es, ap1, tmin

    integer :: it

    tmin = table_ice - 160.
    ap1 = 10. * dim (ta, tmin) + 1.
    ap1 = min (2621., ap1)
    it = ap1
    es = tablew (it) + (ap1 - it) * desw (it)
    wqsat_moist = eps * es * (1. + zvir * qv) / pa

end function wqsat_moist

! =======================================================================
!>@brief The function 'qs1d_m' computes the saturated specific humidity
!! for table iii
! =======================================================================

real function qs1d_m (ta, qv, pa)

    implicit none

    real, intent (in) :: ta, pa, qv

    real :: es, ap1, tmin

    integer :: it

    tmin = table_ice - 160.
    ap1 = 10. * dim (ta, tmin) + 1.
    ap1 = min (2621., ap1)
    it = ap1
    es = table2 (it) + (ap1 - it) * des2 (it)
    qs1d_m = eps * es * (1. + zvir * qv) / pa

end function qs1d_m

! =======================================================================
!>@brief The function 'd_sat' computes the difference in saturation
!! vapor * density * between water and ice
! =======================================================================

real function d_sat (ta, den)

    implicit none

    real, intent (in) :: ta, den

    real :: es_w, es_i, ap1, tmin

    integer :: it

    tmin = table_ice - 160.
    ap1 = 10. * dim (ta, tmin) + 1.
    ap1 = min (2621., ap1)
    it = ap1
    es_w = tablew (it) + (ap1 - it) * desw (it)
    es_i = table2 (it) + (ap1 - it) * des2 (it)
    d_sat = dim (es_w, es_i) / (rvgas * ta * den) ! take positive difference

end function d_sat

! =======================================================================
!>@brief The function 'esw_table' computes the saturated water vapor
!! pressure for table ii
! =======================================================================

real function esw_table (ta)

    implicit none

    real, intent (in) :: ta

    real :: ap1, tmin

    integer :: it

    tmin = table_ice - 160.
    ap1 = 10. * dim (ta, tmin) + 1.
    ap1 = min (2621., ap1)
    it = ap1
    esw_table = tablew (it) + (ap1 - it) * desw (it)

end function esw_table

! =======================================================================
!>@brief The function 'es2_table' computes the saturated water
!! vapor pressure for table iii
! =======================================================================

real function es2_table (ta)

    implicit none

    real, intent (in) :: ta

    real :: ap1, tmin

    integer :: it

    tmin = table_ice - 160.
    ap1 = 10. * dim (ta, tmin) + 1.
    ap1 = min (2621., ap1)
    it = ap1
    es2_table = table2 (it) + (ap1 - it) * des2 (it)

end function es2_table

! =======================================================================
!>\ingroup mod_gfdl_cloud_mp
!! The subroutine 'esw_table1d' computes the saturated water vapor
!! pressure for table ii.
subroutine esw_table1d (ta, es, n)

    implicit none

    integer, intent (in) :: n

    real, intent (in) :: ta (n)

    real, intent (out) :: es (n)

    real :: ap1, tmin

    integer :: i, it

    tmin = table_ice - 160.

    do i = 1, n
        ap1 = 10. * dim (ta (i), tmin) + 1.
        ap1 = min (2621., ap1)
        it = ap1
        es (i) = tablew (it) + (ap1 - it) * desw (it)
    enddo

end subroutine esw_table1d

! =======================================================================
!>\ingroup mod_gfdl_cloud_mp
!! The subroutine 'es3_table1d' computes the saturated water vapor
!! pressure for table iii.
subroutine es2_table1d (ta, es, n)

    implicit none

    integer, intent (in) :: n

    real, intent (in) :: ta (n)

    real, intent (out) :: es (n)

    real :: ap1, tmin

    integer :: i, it

    tmin = table_ice - 160.

    do i = 1, n
        ap1 = 10. * dim (ta (i), tmin) + 1.
        ap1 = min (2621., ap1)
        it = ap1
        es (i) = table2 (it) + (ap1 - it) * des2 (it)
    enddo

end subroutine es2_table1d

! =======================================================================
!>\ingroup mod_gfdl_cloud_mp
!! The subroutine 'es3_table1d' computes the saturated water vapor
!! pressure for table iv.
subroutine es3_table1d (ta, es, n)

    implicit none

    integer, intent (in) :: n

    real, intent (in) :: ta (n)

    real, intent (out) :: es (n)

    real :: ap1, tmin

    integer :: i, it

    tmin = table_ice - 160.

    do i = 1, n
        ap1 = 10. * dim (ta (i), tmin) + 1.
        ap1 = min (2621., ap1)
        it = ap1
        es (i) = table3 (it) + (ap1 - it) * des3 (it)
    enddo

end subroutine es3_table1d

! =======================================================================
!>\ingroup mod_gfdl_cloud_mp
!! saturation water vapor pressure table ii
! 1 - phase table
subroutine qs_tablew (n)

    implicit none

    integer, intent (in) :: n

    real :: delt = 0.1
    real :: tmin, tem, fac0, fac1, fac2

    integer :: i

    tmin = table_ice - 160.

    ! -----------------------------------------------------------------------
    ! compute es over water
    ! -----------------------------------------------------------------------

    do i = 1, n
        tem = tmin + delt * real (i - 1)
        fac0 = (tem - t_ice) / (tem * t_ice)
        fac1 = fac0 * lv0
        fac2 = (dc_vap * log (tem / t_ice) + fac1) / rvgas
        tablew (i) = e00 * exp (fac2)
    enddo

end subroutine qs_tablew

! =======================================================================
!>\ingroup mod_gfdl_cloud_mp
!>@brief saturation water vapor pressure table iii
! 2 - phase table
subroutine qs_table2 (n)

    implicit none

    integer, intent (in) :: n

    real :: delt = 0.1
    real :: tmin, tem0, tem1, fac0, fac1, fac2

    integer :: i, i0, i1

    tmin = table_ice - 160.

    do i = 1, n
        tem0 = tmin + delt * real (i - 1)
        fac0 = (tem0 - t_ice) / (tem0 * t_ice)
        if (i <= 1600) then
            ! -----------------------------------------------------------------------
            ! compute es over ice between - 160 deg c and 0 deg c.
            ! -----------------------------------------------------------------------
            fac1 = fac0 * li2
            fac2 = (d2ice * log (tem0 / t_ice) + fac1) / rvgas
        else
            ! -----------------------------------------------------------------------
            ! compute es over water between 0 deg c and 102 deg c.
            ! -----------------------------------------------------------------------
            fac1 = fac0 * lv0
            fac2 = (dc_vap * log (tem0 / t_ice) + fac1) / rvgas
        endif
        table2 (i) = e00 * exp (fac2)
    enddo

    ! -----------------------------------------------------------------------
    ! smoother around 0 deg c
    ! -----------------------------------------------------------------------

    i0 = 1600
    i1 = 1601
    tem0 = 0.25 * (table2 (i0 - 1) + 2. * table (i0) + table2 (i0 + 1))
    tem1 = 0.25 * (table2 (i1 - 1) + 2. * table (i1) + table2 (i1 + 1))
    table2 (i0) = tem0
    table2 (i1) = tem1

end subroutine qs_table2

! =======================================================================
!>\ingroup mod_gfdl_cloud_mp
!! saturation water vapor pressure table iv
! 2 - phase table with " - 2 c" as the transition point
subroutine qs_table3 (n)

    implicit none

    integer, intent (in) :: n

    real :: delt = 0.1
    real :: esbasw, tbasw, esbasi, tmin, tem, aa, b, c, d, e
    real :: tem0, tem1

    integer :: i, i0, i1

    esbasw = 1013246.0
    tbasw = table_ice + 100.
    esbasi = 6107.1
    tmin = table_ice - 160.

    do i = 1, n
        tem = tmin + delt * real (i - 1)
        ! if (i <= 1600) then
        if (i <= 1580) then ! change to - 2 c
            ! -----------------------------------------------------------------------
            ! compute es over ice between - 160 deg c and 0 deg c.
            ! see smithsonian meteorological tables page 350.
            ! -----------------------------------------------------------------------
            aa = - 9.09718 * (table_ice / tem - 1.)
            b = - 3.56654 * alog10 (table_ice / tem)
            c = 0.876793 * (1. - tem / table_ice)
            e = alog10 (esbasi)
            table3 (i) = 0.1 * 10 ** (aa + b + c + e)
        else
            ! -----------------------------------------------------------------------
            ! compute es over water between - 2 deg c and 102 deg c.
            ! see smithsonian meteorological tables page 350.
            ! -----------------------------------------------------------------------
            aa = - 7.90298 * (tbasw / tem - 1.)
            b = 5.02808 * alog10 (tbasw / tem)
            c = - 1.3816e-7 * (10 ** ((1. - tem / tbasw) * 11.344) - 1.)
            d = 8.1328e-3 * (10 ** ((tbasw / tem - 1.) * (- 3.49149)) - 1.)
            e = alog10 (esbasw)
            table3 (i) = 0.1 * 10 ** (aa + b + c + d + e)
        endif
    enddo

    ! -----------------------------------------------------------------------
    ! smoother around - 2 deg c
    ! -----------------------------------------------------------------------

    i0 = 1580
    i1 = 1581
    tem0 = 0.25 * (table3 (i0 - 1) + 2. * table (i0) + table3 (i0 + 1))
    tem1 = 0.25 * (table3 (i1 - 1) + 2. * table (i1) + table3 (i1 + 1))
    table3 (i0) = tem0
    table3 (i1) = tem1

end subroutine qs_table3

! =======================================================================
! compute the saturated specific humidity for table
! note: this routine is based on "moist" mixing ratio
!>\ingroup mod_gfdl_cloud_mp
!! The function 'qs_blend' computes the saturated specific humidity
!! with a blend of water and ice depending on the temperature.
real function qs_blend (t, p, q)

    implicit none

    real, intent (in) :: t, p, q

    real :: es, ap1, tmin

    integer :: it

    tmin = table_ice - 160.
    ap1 = 10. * dim (t, tmin) + 1.
    ap1 = min (2621., ap1)
    it = ap1
    es = table (it) + (ap1 - it) * des (it)
    qs_blend = eps * es * (1. + zvir * q) / p

end function qs_blend

! =======================================================================
!>\ingroup mod_gfdl_cloud_mp
!! saturation water vapor pressure table i
! 3 - phase table
subroutine qs_table (n)

    implicit none

    integer, intent (in) :: n

    real :: delt = 0.1
    real :: tmin, tem, esh20
    real :: wice, wh2o, fac0, fac1, fac2
    real :: esupc (200)

    integer :: i

    tmin = table_ice - 160.

    ! -----------------------------------------------------------------------
    ! compute es over ice between - 160 deg c and 0 deg c.
    ! -----------------------------------------------------------------------

    do i = 1, 1600
        tem = tmin + delt * real (i - 1)
        fac0 = (tem - t_ice) / (tem * t_ice)
        fac1 = fac0 * li2
        fac2 = (d2ice * log (tem / t_ice) + fac1) / rvgas
        table (i) = e00 * exp (fac2)
    enddo

    ! -----------------------------------------------------------------------
    ! compute es over water between - 20 deg c and 102 deg c.
    ! -----------------------------------------------------------------------

    do i = 1, 1221
        tem = 253.16 + delt * real (i - 1)
        fac0 = (tem - t_ice) / (tem * t_ice)
        fac1 = fac0 * lv0
        fac2 = (dc_vap * log (tem / t_ice) + fac1) / rvgas
        esh20 = e00 * exp (fac2)
        if (i <= 200) then
            esupc (i) = esh20
        else
            table (i + 1400) = esh20
        endif
    enddo

    ! -----------------------------------------------------------------------
    ! derive blended es over ice and supercooled water between - 20 deg c and 0 deg c
    ! -----------------------------------------------------------------------

    do i = 1, 200
        tem = 253.16 + delt * real (i - 1)
        wice = 0.05 * (table_ice - tem)
        wh2o = 0.05 * (tem - 253.16)
        table (i + 1400) = wice * table (i + 1400) + wh2o * esupc (i)
    enddo

end subroutine qs_table

! =======================================================================
! compute the saturated specific humidity and the gradient of saturated specific humidity
! input t in deg k, p in pa; p = rho rdry tv, moist pressure
!>\ingroup mod_gfdl_cloud_mp
!! The function 'qsmith' computes the saturated specific humidity
!! with a blend of water and ice depending on the temperature in 3D.
!@details It als oincludes the option for computing des/dT.
subroutine qsmith (im, km, ks, t, p, q, qs, dqdt)

    implicit none

    integer, intent (in) :: im, km, ks

    real, intent (in), dimension (im, km) :: t, p, q

    real, intent (out), dimension (im, km) :: qs

    real, intent (out), dimension (im, km), optional :: dqdt

    real :: eps10, ap1, tmin

    real, dimension (im, km) :: es

    integer :: i, k, it

    tmin = table_ice - 160.
    eps10 = 10. * eps

    if (.not. tables_are_initialized) then
        call qsmith_init
    endif

    do k = ks, km
        do i = 1, im
            ap1 = 10. * dim (t (i, k), tmin) + 1.
            ap1 = min (2621., ap1)
            it = ap1
            es (i, k) = table (it) + (ap1 - it) * des (it)
            qs (i, k) = eps * es (i, k) * (1. + zvir * q (i, k)) / p (i, k)
        enddo
    enddo

    if (present (dqdt)) then
        do k = ks, km
            do i = 1, im
                ap1 = 10. * dim (t (i, k), tmin) + 1.
                ap1 = min (2621., ap1) - 0.5
                it = ap1
                dqdt (i, k) = eps10 * (des (it) + (ap1 - it) * (des (it + 1) - des (it))) * (1. + zvir * q (i, k)) / p (i, k)
            enddo
        enddo
    endif

end subroutine qsmith

! =======================================================================
!>\ingroup mod_gfdl_cloud_mp
!>@brief The subroutine 'neg_adj' fixes negative water species.
!>@details This is designed for 6-class micro-physics schemes.
subroutine neg_adj (ktop, kbot, pt, dp, qv, ql, qr, qi, qs, qg)

    implicit none

    integer, intent (in) :: ktop, kbot

    real, intent (in), dimension (ktop:kbot) :: dp

    real, intent (inout), dimension (ktop:kbot) :: pt, qv, ql, qr, qi, qs, qg

    real, dimension (ktop:kbot) :: lcpk, icpk

    real :: dq, cvm

    integer :: k

    ! -----------------------------------------------------------------------
    ! define heat capacity and latent heat coefficient
    ! -----------------------------------------------------------------------

    do k = ktop, kbot
        cvm = c_air + qv (k) * c_vap + (qr (k) + ql (k)) * c_liq + (qi (k) + qs (k) + qg (k)) * c_ice
        lcpk (k) = (lv00 + d0_vap * pt (k)) / cvm
        icpk (k) = (li00 + dc_ice * pt (k)) / cvm
    enddo

    do k = ktop, kbot

        ! -----------------------------------------------------------------------
        ! ice phase:
        ! -----------------------------------------------------------------------

        ! if cloud ice < 0, borrow from snow
        if (qi (k) < 0.) then
            qs (k) = qs (k) + qi (k)
            qi (k) = 0.
        endif
        ! if snow < 0, borrow from graupel
        if (qs (k) < 0.) then
            qg (k) = qg (k) + qs (k)
            qs (k) = 0.
        endif
        ! if graupel < 0, borrow from rain
        if (qg (k) < 0.) then
            qr (k) = qr (k) + qg (k)
            pt (k) = pt (k) - qg (k) * icpk (k) ! heating
            qg (k) = 0.
        endif

        ! -----------------------------------------------------------------------
        ! liquid phase:
        ! -----------------------------------------------------------------------

        ! if rain < 0, borrow from cloud water
        if (qr (k) < 0.) then
            ql (k) = ql (k) + qr (k)
            qr (k) = 0.
        endif
        ! if cloud water < 0, borrow from water vapor
        if (ql (k) < 0.) then
            qv (k) = qv (k) + ql (k)
            pt (k) = pt (k) - ql (k) * lcpk (k) ! heating
            ql (k) = 0.
        endif

    enddo

    ! -----------------------------------------------------------------------
    ! fix water vapor; borrow from below
    ! -----------------------------------------------------------------------

    do k = ktop, kbot - 1
        if (qv (k) < 0.) then
            qv (k + 1) = qv (k + 1) + qv (k) * dp (k) / dp (k + 1)
            qv (k) = 0.
        endif
    enddo

    ! -----------------------------------------------------------------------
    ! bottom layer; borrow from above
    ! -----------------------------------------------------------------------

    if (qv (kbot) < 0. .and. qv (kbot - 1) > 0.) then
        dq = min (- qv (kbot) * dp (kbot), qv (kbot - 1) * dp (kbot - 1))
        qv (kbot - 1) = qv (kbot - 1) - dq / dp (kbot - 1)
        qv (kbot) = qv (kbot) + dq / dp (kbot)
    endif

end subroutine neg_adj

! =======================================================================
! compute global sum
!>@brief quick local sum algorithm
! =======================================================================

!real function g_sum (p, ifirst, ilast, jfirst, jlast, area, mode)
!
! use mpp_mod, only: mpp_sum
!
! implicit none
!
! integer, intent (in) :: ifirst, ilast, jfirst, jlast
! integer, intent (in) :: mode ! if == 1 divided by area
!
! real, intent (in), dimension (ifirst:ilast, jfirst:jlast) :: p, area
!
! integer :: i, j
!
! real :: gsum
!
! if (global_area < 0.) then
! global_area = 0.
! do j = jfirst, jlast
! do i = ifirst, ilast
! global_area = global_area + area (i, j)
! enddo
! enddo
! call mpp_sum (global_area)
! endif
!
! gsum = 0.
! do j = jfirst, jlast
! do i = ifirst, ilast
! gsum = gsum + p (i, j) * area (i, j)
! enddo
! enddo
! call mpp_sum (gsum)
!
! if (mode == 1) then
! g_sum = gsum / global_area
! else
! g_sum = gsum
! endif
!
!end function g_sum

! ==========================================================================
!>\ingroup mod_gfdl_cloud_mp
!! The subroutine 'interpolate_z' interpolates to a prescribed height.
subroutine interpolate_z (is, ie, js, je, km, zl, hgt, a3, a2)

    implicit none

    integer, intent (in) :: is, ie, js, je, km

    real, intent (in), dimension (is:ie, js:je, km) :: a3

    real, intent (in), dimension (is:ie, js:je, km + 1) :: hgt ! hgt (k) > hgt (k + 1)

    real, intent (in) :: zl

    real, intent (out), dimension (is:ie, js:je) :: a2

    real, dimension (km) :: zm !< middle layer height

    integer :: i, j, k

    !$omp parallel do default (none) shared (is, ie, js, je, km, hgt, zl, a2, a3) private (zm)

    do j = js, je
        do i = is, ie
            do k = 1, km
                zm (k) = 0.5 * (hgt (i, j, k) + hgt (i, j, k + 1))
            enddo
            if (zl >= zm (1)) then
                a2 (i, j) = a3 (i, j, 1)
            elseif (zl <= zm (km)) then
                a2 (i, j) = a3 (i, j, km)
            else
                do k = 1, km - 1
                    if (zl <= zm (k) .and. zl >= zm (k + 1)) then
                        a2 (i, j) = a3 (i, j, k) + (a3 (i, j, k + 1) - a3 (i, j, k)) * (zm (k) - zl) / (zm (k) - zm (k + 1))
                        exit
                    endif
                enddo
            endif
        enddo
    enddo

end subroutine interpolate_z

! =======================================================================
!> \ingroup mod_gfdl_cloud_mp
!! The subroutine 'cloud_diagnosis' diagnoses the radius of cloud
!! species.
!>\author Linjiong Zhoum, Shian-Jiann Lin
! =======================================================================
subroutine cloud_diagnosis (is, ie, ks, ke, den, delp, lsm, qmw, qmi, qmr, qms, qmg, t, &
        rew, rei, rer, res, reg)

    implicit none

    integer, intent (in) :: is, ie, ks, ke
    integer, intent (in), dimension (is:ie) :: lsm ! land sea mask, 0: ocean, 1: land, 2: sea ice

    real, intent (in), dimension (is:ie, ks:ke) :: den, delp, t
    real, intent (in), dimension (is:ie, ks:ke) :: qmw, qmi, qmr, qms, qmg !< units: kg / kg

    real, intent (out), dimension (is:ie, ks:ke) :: rew, rei, rer, res, reg !< units: micron

    real, dimension (is:ie, ks:ke) :: qcw, qci, qcr, qcs, qcg !< units: g / m^2
    
    integer :: i, k

    real :: lambdar, lambdas, lambdag
    real :: dpg, rei_fac, mask, ccn, bw
    real, parameter :: rho_0 = 50.e-3

    real :: rhow = 1.0e3, rhor = 1.0e3, rhos = 1.0e2, rhog = 4.0e2
    real :: n0r = 8.0e6, n0s = 3.0e6, n0g = 4.0e6
    real :: alphar = 0.8, alphas = 0.25, alphag = 0.5
    real :: gammar = 17.837789, gammas = 8.2850630, gammag = 11.631769
    real :: qmin = 1.0e-12, beta = 1.22, qmin1 = 9.e-6

    do k = ks, ke
        do i = is, ie
            
            dpg = abs (delp (i, k)) / grav
            mask = min (max (real(lsm (i)), 0.0), 2.0)

            ! -----------------------------------------------------------------------
            ! cloud water (Martin et al., 1994)
            ! -----------------------------------------------------------------------

            ccn = 0.80 * (- 1.15e-3 * (ccn_o ** 2) + 0.963 * ccn_o + 5.30) * abs (mask - 1.0) + &
                0.67 * (- 2.10e-4 * (ccn_l ** 2) + 0.568 * ccn_l - 27.9) * (1.0 - abs (mask - 1.0))

            if (qmw (i, k) .gt. qmin) then
                qcw (i, k) = dpg * qmw (i, k) * 1.0e3
                rew (i, k) = exp (1.0 / 3.0 * log ((3.0 * den (i, k) * qmw (i, k)) / (4.0 * pi * rhow * ccn))) * 1.0e4
                rew (i, k) = max (rewmin, min (rewmax, rew (i, k)))
            else
                qcw (i, k) = 0.0
                rew (i, k) = rewmin
            endif
            
            if (reiflag .eq. 1) then

            ! -----------------------------------------------------------------------
                ! cloud ice (Heymsfield and Mcfarquhar, 1996)
            ! -----------------------------------------------------------------------

                if (qmi (i, k) .gt. qmin1) then
                    qci (i, k) = dpg * qmi (i, k) * 1.0e3
                    rei_fac = log (1.0e3 * qmi (i, k) * den (i, k))
                    if (t (i, k) - tice .lt. - 50) then
                        rei (i, k) = beta / 9.917 * exp (0.109 * rei_fac) * 1.0e3
                    elseif (t (i, k) - tice .lt. - 40) then
                        rei (i, k) = beta / 9.337 * exp (0.080 * rei_fac) * 1.0e3
                    elseif (t (i, k) - tice .lt. - 30) then
                        rei (i, k) = beta / 9.208 * exp (0.055 * rei_fac) * 1.0e3
                else
                        rei (i, k) = beta / 9.387 * exp (0.031 * rei_fac) * 1.0e3
                endif
                    rei (i, k) = max (reimin, min (reimax, rei (i, k)))
            else
                    qci (i, k) = 0.0
                    rei (i, k) = reimin
            endif

            endif
            
            if (reiflag .eq. 2) then

            ! -----------------------------------------------------------------------
                ! cloud ice (Wyser, 1998)
            ! -----------------------------------------------------------------------

                if (qmi (i, k) .gt. qmin1) then
                    qci (i, k) = dpg * qmi (i, k) * 1.0e3
                    bw = - 2. + 1.e-3 * log10 (den (i, k) * qmi (i, k) / rho_0) * max (0.0, tice - t (i, k)) ** 1.5
                    rei (i, k) = 377.4 + bw * (203.3 + bw * (37.91 + 2.3696 * bw))
                    rei (i, k) = max (reimin, min (reimax, rei (i, k)))
            else
                    qci (i, k) = 0.0
                    rei (i, k) = reimin
                endif

            endif

            ! -----------------------------------------------------------------------
            ! rain (Lin et al., 1983)
            ! -----------------------------------------------------------------------

            if (qmr (i, k) .gt. qmin) then
                qcr (i, k) = dpg * qmr (i, k) * 1.0e3
                lambdar = exp (0.25 * log (pi * rhor * n0r / qmr (i, k) / den (i, k)))
                rer (i, k) = 0.5 * exp (log (gammar / 6) / alphar) / lambdar * 1.0e6
                rer (i, k) = max (rermin, min (rermax, rer (i, k)))
            else
                qcr (i, k) = 0.0
                rer (i, k) = rermin
            endif

            ! -----------------------------------------------------------------------
            ! snow (Lin et al., 1983)
            ! -----------------------------------------------------------------------

            if (qms (i, k) .gt. qmin1) then
                qcs (i, k) = dpg * qms (i, k) * 1.0e3
                lambdas = exp (0.25 * log (pi * rhos * n0s / qms (i, k) / den (i, k)))
                res (i, k) = 0.5 * exp (log (gammas / 6) / alphas) / lambdas * 1.0e6
                res (i, k) = max (resmin, min (resmax, res (i, k)))
            else
                qcs (i, k) = 0.0
                res (i, k) = resmin
            endif
            
            ! -----------------------------------------------------------------------
            ! graupel (Lin et al., 1983)
            ! -----------------------------------------------------------------------
            
            if (qmg (i, k) .gt. qmin) then
                qcg (i, k) = dpg * qmg (i, k) * 1.0e3
                lambdag = exp (0.25 * log (pi * rhog * n0g / qmg (i, k) / den (i, k)))
                reg (i, k) = 0.5 * exp (log (gammag / 6) / alphag) / lambdag * 1.0e6
                reg (i, k) = max (regmin, min (regmax, reg (i, k)))
            else
                qcg (i, k) = 0.0
                reg (i, k) = regmin
            endif

        enddo
    enddo

end subroutine cloud_diagnosis

!+---+-----------------------------------------------------------------+
!>\ingroup mod_gfdl_cloud_mp
!! This subroutine calculates radar reflectivity.
      subroutine refl10cm_gfdl (qv1d, qr1d, qs1d, qg1d,                 &
                       t1d, p1d, dBZ, kts, kte, ii, jj, melti)

      IMPLICIT NONE

!..Sub arguments
      INTEGER, INTENT(IN):: kts, kte, ii,jj
      REAL, DIMENSION(kts:kte), INTENT(IN)::                            &
                      qv1d, qr1d, qs1d, qg1d, t1d, p1d
      REAL, DIMENSION(kts:kte), INTENT(INOUT):: dBZ

!..Local variables
      REAL, DIMENSION(kts:kte):: temp, pres, qv, rho
      REAL, DIMENSION(kts:kte):: rr, rs, rg
!      REAL:: temp_C

      DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilams, ilamg
      DOUBLE PRECISION, DIMENSION(kts:kte):: N0_r, N0_s, N0_g
      DOUBLE PRECISION:: lamr, lams, lamg
      LOGICAL, DIMENSION(kts:kte):: L_qr, L_qs, L_qg

      REAL, DIMENSION(kts:kte):: ze_rain, ze_snow, ze_graupel
      DOUBLE PRECISION:: fmelt_s, fmelt_g

      INTEGER:: i, k, k_0, kbot, n
      LOGICAL, INTENT(IN):: melti
      DOUBLE PRECISION:: cback, x, eta, f_d
!+---+

      do k = kts, kte
         dBZ(k) = -35.0
      enddo

!+---+-----------------------------------------------------------------+
!..Put column of data into local arrays.
!+---+-----------------------------------------------------------------+
      do k = kts, kte
         temp(k) = t1d(k)
!         temp_C = min(-0.001, temp(K)-273.15)
         qv(k) = MAX(1.E-10, qv1d(k))
         pres(k) = p1d(k)
         rho(k) = 0.622*pres(k)/(rdgas*temp(k)*(qv(k)+0.622))

         if (qr1d(k) .gt. 1.E-9) then
            rr(k) = qr1d(k)*rho(k)
            N0_r(k) = n0r
            lamr = (xam_r*xcrg(3)*N0_r(k)/rr(k))**(1./xcre(1))
            ilamr(k) = 1./lamr
            L_qr(k) = .true.
         else
            rr(k) = 1.E-12
            L_qr(k) = .false.
         endif

         if (qs1d(k) .gt. 1.E-9) then
            rs(k) = qs1d(k)*rho(k)
            N0_s(k) = n0s
            lams = (xam_s*xcsg(3)*N0_s(k)/rs(k))**(1./xcse(1))
            ilams(k) = 1./lams
            L_qs(k) = .true.
         else
            rs(k) = 1.E-12
            L_qs(k) = .false.
         endif

         if (qg1d(k) .gt. 1.E-9) then
            rg(k) = qg1d(k)*rho(k)
            N0_g(k) = n0g
            lamg = (xam_g*xcgg(3)*N0_g(k)/rg(k))**(1./xcge(1))
            ilamg(k) = 1./lamg
            L_qg(k) = .true.
         else
            rg(k) = 1.E-12
            L_qg(k) = .false.
         endif
      enddo

!+---+-----------------------------------------------------------------+
!..Locate K-level of start of melting (k_0 is level above).
!+---+-----------------------------------------------------------------+
      k_0 = kts
      K_LOOP:do k = kte-1, kts, -1
         if ( melti .and. (temp(k).gt.273.15) .and. L_qr(k)             &
              .and. (L_qs(k+1).or.L_qg(k+1)) ) then
            k_0 = MAX(k+1, k_0)
            EXIT K_LOOP
         endif
      enddo K_LOOP
!+---+-----------------------------------------------------------------+
!..Assume Rayleigh approximation at 10 cm wavelength. Rain (all temps)
!.. and non-water-coated snow and graupel when below freezing are
!.. simple. Integrations of m(D)*m(D)*N(D)*dD.
!+---+-----------------------------------------------------------------+
      do k = kts, kte
         ze_rain(k) = 1.e-22
         ze_snow(k) = 1.e-22
         ze_graupel(k) = 1.e-22
         if (L_qr(k)) ze_rain(k) = N0_r(k)*xcrg(4)*ilamr(k)**xcre(4)
         if (L_qs(k)) ze_snow(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI)     &
                                 * (xam_s/900.0)*(xam_s/900.0)          &
                                 * N0_s(k)*xcsg(4)*ilams(k)**xcse(4)
         if (L_qg(k)) ze_graupel(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI)  &
                                    * (xam_g/900.0)*(xam_g/900.0)       &
                                    * N0_g(k)*xcgg(4)*ilamg(k)**xcge(4)
      enddo


!+---+-----------------------------------------------------------------+
!..Special case of melting ice (snow/graupel) particles.  Assume the
!.. ice is surrounded by the liquid water.  Fraction of meltwater is
!.. extremely simple based on amount found above the melting level.
!.. Uses code from Uli Blahak (rayleigh_soak_wetgraupel and supporting
!.. routines).
!+---+-----------------------------------------------------------------+

      if (melti .and. k_0.ge.kts+1) then
       do k = k_0-1, kts, -1

!..Reflectivity contributed by melting snow
          if (L_qs(k) .and. L_qs(k_0) ) then
           fmelt_s = MAX(0.005d0, MIN(1.0d0-rs(k)/rs(k_0), 0.99d0))
           eta = 0.d0
           lams = 1./ilams(k)
           do n = 1, nrbins
              x = xam_s * xxDs(n)**xbm_s
              call rayleigh_soak_wetgraupel (x,DBLE(xocms),DBLE(xobms), &
                    fmelt_s, melt_outside_s, m_w_0, m_i_0, lamda_radar, &
                    CBACK, mixingrulestring_s, matrixstring_s,          &
                    inclusionstring_s, hoststring_s,                    &
                    hostmatrixstring_s, hostinclusionstring_s)
              f_d = N0_s(k)*xxDs(n)**xmu_s * DEXP(-lams*xxDs(n))
              eta = eta + f_d * CBACK * simpson(n) * xdts(n)
           enddo
           ze_snow(k) = SNGL(lamda4 / (pi5 * K_w) * eta)
          endif


!..Reflectivity contributed by melting graupel

          if (L_qg(k) .and. L_qg(k_0) ) then
           fmelt_g = MAX(0.005d0, MIN(1.0d0-rg(k)/rg(k_0), 0.99d0))
           eta = 0.d0
           lamg = 1./ilamg(k)
           do n = 1, nrbins
              x = xam_g * xxDg(n)**xbm_g
              call rayleigh_soak_wetgraupel (x,DBLE(xocmg),DBLE(xobmg), &
                    fmelt_g, melt_outside_g, m_w_0, m_i_0, lamda_radar, &
                    CBACK, mixingrulestring_g, matrixstring_g,          &
                    inclusionstring_g, hoststring_g,                    &
                    hostmatrixstring_g, hostinclusionstring_g)
              f_d = N0_g(k)*xxDg(n)**xmu_g * DEXP(-lamg*xxDg(n))
              eta = eta + f_d * CBACK * simpson(n) * xdtg(n)
           enddo
           ze_graupel(k) = SNGL(lamda4 / (pi5 * K_w) * eta)
          endif

       enddo
      endif

      do k = kte, kts, -1
         dBZ(k) = 10.*log10((ze_rain(k)+ze_snow(k)+ze_graupel(k))*1.d18)
      enddo


      end subroutine refl10cm_gfdl
!+---+-----------------------------------------------------------------+

end module gfdl_cloud_microphys_mod
